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1.  Introduction 


This  second  semi-annual  progress  report  contains  a  summary  of  work  accomplished  on 
O.N.R.  contract  number  N00014-86-K.-0370,  Delay-Doppler  Radar-Imaging,  during  the  period 
from  1  December  1986  to  30  May  1987. 

The  goal  of  this  project  is  to  formulate  and  investigate  new  approaches  for  forming  images 
of  radar  targets  from  spotlight- mode,  delay-doppler  measurements.  These  measurements  could 
be  acquired  with  a  high- resolution  radar-imaging  system  operating  with  an  optical  or  radio 
frequency  carrier.  Two  approaches  are  under  study.  The  first  is  motivated  by  an 
image- reconstruction  algorithm  used  in  radionuclide  imaging  called  the  "confidence-weighted” 
algorithm.  The  second  is  one  based  on  more  fundamental  principles  which  starts  with  a 
mathematical  model  that  accurately  describes  the  physics  of  an  imaging  radar-system  and  then 
uses  statistical-estimation  theory  with  this  model  to  derive  processing  algorithms. 

Spotlight- mode  high-resolution  radar- imaging  relies  upon  the  relative  motion  between  the 
transmitter,  target,  and  receiver.  The  target  is  illuminated  by  a  series  of  transmitted  pulses. 
The  return  for  each  pulse  is  a  superposition  of  reflections  from  various  locations  on  the  target, 
with  each  location  affecting  the  pulse  by  introducing  both  a  delay  and  a  doppler  shift.  The 
returns  are  processed  to  produce  an  image  of  the  target.  ■4 

The  common  approach  is  to  use  the  same  transmitted  pulse  for  each  illumination  of  the 
target.  The  returns  are  processed  using  a  two-dimensional  Fourier  transform  to  produce  the 
target’s  image  in  delay /doppler  or  range/cross -range  coordinates  [1,2].  Our  goal  is  to  compare 
images  produced  in  this  standard  way  with  those  produced  using  the  alternative  approaches  we 
are  developing. 

Bernfeld  [3]  appears  to  be  the  first  to  introduce  the  idea  for  radar  imaging  of  modifying 
the  pulse  shape  on  successive  illuminations  of  the  target.  We  are  using  Bernfeld’s  idea  in  the 
confidence- weighted  approach.  With  this,  the  FM  chirp  rate  of  each  pulse  is  varied  so  that  the 
angles  made  by  the  ambiguity  functions  in  the  delay /doppler  plane  are  caused  to  vary  over  the 
full  range  of  angles  between  0*  and  180*.  Use  is  then  made  of  the  fact  that,  on  the  average,  the 
output  of  a  receiver  consisting  of  filters  matched  to  doppler  shifted  versions  of  a  transmitted 
pulse  is  the  two-dimensional  convolution  of  ambiguity  function  of  the  pulse  with  the  scattering 


function  of  the  target  (4],  an  output  we  call  the  delay/doppler  power  function.  Given  the 
delay/doppler  power  functions  for  each  illumination,  the  target’s  ambiguity  function  can  be 
determined  using  the  confidence-weighted  algorithm  [5]. 

2.  Summary  of  Work  Accomplished 

During  this  reporting  period,  Mr.  Robert  C.  Lewis  completed  the  implementation  of 
conventional,  two-dimensional  Fourier- transform  processing  for  use  in  comparison  studies  with 
the  new  algorithms  we  are  developing.  This  work  is  documented  in  his  M.S.  thesis  [6],  which  is 
included  here  in  Appendix  4.1.  A  computer  simulation  was  developed  of  inverse  synthet¬ 
ic-aperture  radar  processing  using  a  stepped -frequency  transmitted  waveform.  This  simulation 
was  validated  by  the  generation  of  several  images  from  target  scattering  functions.  A 
mathematical  model  for  a  slowly  fluctuating  point  target  in  the  presence  of  additive  noise  was 
applied  to  simulating  the  received  signal  from  an  extended  diffuse  target.  This  model  was 
implemented  in  the  computer  simulation.  The  effects  on  conventional  radar- images  of  the 
quantitatively  specified  noise  are  given.  Images  of  simple  radar  targets  were  produced  for  a 
random  reflectivity  process  having  varying  degrees  of  coherence  time.  The  results  show  that 
random  variations  in  the  reflectivity  have  a  significant  effect  on  the  quality  of  conventional 
ISAR  images.  The  target  is  unrecognizable  unless  the  reflectivity  has  a  long  coherence  time 
compared  to  the  duration  of  the  radar  illumination. 

We  are  currently  performing  a  similar  simulation  in  which  the  confidence-weighted 
approach  is  used  to  produce  the  images.  The  results  will  be  presented  in  the  next  progress 
report.  A  preliminary  observation  is  that  the  images  produced  this  way  appear  to  be  less 
degraded  by  random  temporal  variations  in  the  reflectivity. 

The  use  of  the  two-dimensional  Fourier  transform  to  produce  target  images  is  based  on  a 
deterministic  model  for  the  radar -reflection  data.  However,  commonly  accepted  models  for 
returns  from  diffuse  targets  at  optical  and  radio  frequencies  is  not  deterministic.  We  anticipate 
that  improved  target  images  will  be  obtained  by  taking  the  random  nature  of  the  radar  return 
into  account  during  the  processing.  For  this  reason,  we  have  initiated  work  in  formulating  the 
high-resolution  radar-imaging  problem  as  a  problem  in  statistical  estimation.  Our  first  results 
have  been  prepared  as  a  report  [7],  which  is  included  here  in  Appendix  4.2.  This  report  is 
presently  being  revised  for  publication  in  the  IEEE  Trans,  on  Information  Theory,  with  the 


main  objective  in  the  revision  being  to  reduce  the  dimensionality  of  the  matrices  required  to 
form  the  image.  The  main  result  of  this  report  is  an  iterative  algorithm  for  producing  target 
images. 
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INVERSE  SYNTHETIC -APERTURE  RADAR  IMAGING 


1.  INTRODUCTION 


Inverse  synthetic-aperture  radar  (ISAR)  imaging  is  a  technique  for 
making  two-dimensional  delay-doppler  images  of  flying  aircraft  at  long 
distances.  The  resulting  image  is  a  spatial  distribution  of  the 
reflected  signal  power  from  the  surface  of  the  target  at  sufficient 
resolution  to  distinguish  features  for  recognition.  These  radar  images 
are  often  judged  similar  enough  to  visual  images  to  enable  feature 
recognition  and  target  identification. 

There  are  many  other  applications  in  which  the  ISAR  imaging 
techniques  are  used.  Sensors  using  other  frequencies  and  signal  media 
are  used;  for  example,  underwater  sound  (SONAR)  is  used  in  submarines  to 
make  images  of  other  ships.  Also,  the  idea  of  forming  an  image  of  an 
unknown  target  for  recognition  purposes  is  not  the  only  objective. 
Radar  images  are  made  in  laboratories  to  resolve  the  radar  cross  section 
of  features  on  an  aircraft  or  other  radar  reflector. 

This  thesis  addresses  the  application  of  airborne  imaging  of  a 
target  aircraft.  It  documents  the  results  of  a  project  during  which  a 
computer  program  was  developed  to  generate  an  ISAR  image  from  a  speci¬ 
fied  target  scattering  function.  The  results  of  this  development  are 
destonstrated  by  exasiples  of  images  which  were  generated.  In  addition,  a 
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noise  model  was  developed  for  radar  data  and  was  implemented  in  the 
simulation  to  provide  a  quantitative  analysis  tool  useful  for  observing 
the  effects  of  random  reflectivity  and  additive  channel  noise  on  radar 
images.  Resulting  images  are  presented. 

1.1  BACKGROUND 

The  research  for  this  thesis  was  conducted  as  a  part  of  a  radar 
imaging  project  in  the  newly  formed  Electronic  Systems  and  Signals 
Research  Laboratory  of  the  Electrical  Engineering  department  at 
Washington  University.  It  contributes  to  the  laboratory's  objective  of 
developing  an  ISAR  processing  capability  to  be  used  as  an  aid  in  study* 
ing  the  performance  of  conventional  ISAR  imaging  compared  to  new  imaging 
algorithms  currently  in  development,  see  (1)  .  The  new  algorithms  are 
based  on  different  approaches  to  the  radar  imaging  problem,  and  their 
relative  performance  to  ISAR  is  a  major  laboratory  interest. 

1.2  RADAR  IMAGING 

The  signal  transmitted  by  a  radar  can  be  designed  so  that  received 
echo  signals  can  be  measured  to  yield  information  about  the  appearance 
of  the  target.  The  measurements  are  usually  displayed  as  a  two-dimen¬ 
sional  image,  like  a  photograph.  A  well-known  experience  which  is  an 
analogy  to  this  operation  is  using  a  flashlight  in  a  dark  room  to 
illuminate  an  object  within  the  light  beam.  The  flashlight  is  analogous 


*  The  numbers  in  parentheses  in  the  text  indicate  references  in  the 
Bibliography. 
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to  the  radar  illuminating  the  target.  A  radar  illuminates  a  target  with 
microwave  energy  instead  of  light.  The  eye  can  resolve  the  illuminated 
surface  of  an  object  quite  finely  so  that  detailed  features  are  seen. 
Likewise  the  radar  processing  must  finely  resolve  the  surface  of  the 
illuminated  target.  An  image  shows  the  variation  of  reflected  power 
across  the  surface  of  the  object  to  the  spatial  resolution  limit.  The 
radar's  target  also  resides  within  the  beam  width  of  illumination.  One 
important  difference  of  radar  imaging,  compared  with  visual  imaging,  is 
that  the  target  must  be  rotating. 

The  radar-signal  processing  mist  make  use  of  two  received  signal 
parameters  in  order  to  produce  an  image.  The  paraamters  must  be 
directly  related  to  separate  orthogonal  components  of  the  two-dimen¬ 
sional  position  across  the  object  as  projected  on  the  image  plane.  In 
radar  imaging,  the  received  signal  parameters  used  are  round-trip  delay 
and  doppler  frequency  shift.  Any  point  on  the  surface  of  the  target  can 
be  located  in  a  delay-doppler  two-dimensional  image  plane.  To  under¬ 
stand  this,  refer  to  the  geometry  of  the  location  of  the  radar  set 
relative  to  the  target,  as  shorn  in  Figure  1.1.  A  narrow  pulse  trans¬ 
mitted  from  the  radar  and  reflected  from  a  point  target  will  be  received 
after  a  round-trip  delay  time.  The  delay  locates  the  point  in 
down-range  position.  A  large  target  is  a  rigid  body  rotating  about  an 
axis  at  its  geometric  center,  and  points  distributed  on  its  surface  are 
moving  at  different  velocities,  depending  on  each  one's  distance  from 
the  axis.  For  a  group  of  points  with  the  same  delay  location,  contained 
in  a  slice  of  the  target  at  a  down-range  position,  their  cross-range 
position  corresponds  to  distance  from  the  axis  of  rotation,  and  there¬ 
fore  velocity.  A  target  point's  velocity  causes  a  doppler  frequency 
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shift  of  the  reflected  radar  signal.  By  Manuring  dopplar  frequency,  a 
point  target  can  thus  be  located  in  cross-range  by  using  doppler.  The 
delay-doppler  plane  indicated  in  Figure  1.1  corresponds  directly  to  the 
x-y  position  plane  shown. 


1.3  ISA*  IMAGING 

The  relative  rotation  between  the  radar  and  the  target  may  come 
from  either  the  radar  being  fixed  while  the  target  rotates,  or  equi¬ 
valently,  the  radar  rotates  about  the  fixed  target.  The  latter  situ¬ 
ation  is  used  in  synthetic  aperture  radar  (SAR)  imaging  where  a  ground 
map  is  made  as  an  aircraft  flies  by.  The  descriptive  term  "synthetic 
aperture"  refers  to  the  equivalent  antenna  aperture  derived  from 


integrating  received  signals  over  a  long  flight  path.  A  longer  antenna 
results  in  finer  resolution  of  the  target,  thus  a  long  flight  path  while 
imaging  is  desirable.  Ground  sapping  by  SAR  techniques  preceded  the  use 
of  ISAR,  or  inverse  SAR,  for  aircraft- to-aircraft  iaaging.  The  descrip¬ 
tion  "inverse"  refers  to  using  a  fixed  radar  and  a  rotating  target, 
which  is  close  to  the  actual  situation  in  air-to-air  iaaging  when  using 
the  radar's  aircraft  as  a  position  reference.  The  cooputer  processing 
used  in  the  above  two  cases  of  iaaging  are  nearly  the  same. 

The  ISAR  imaging  geoaetry  is  shown  in  Figure  1.1.  The  axis  of 
rotation  is  normal  to  the  line  between  the  radar  and  the  center  of  the 
target.  The  resulting  image  shows  the  point  of  view  of  looking  along 
the  axis  of  rotation,  with  illuaination  coming  from  the  side. 

A  large  target's  body  rotation  axis  may  have  arbitrary  angle  of 
intersection  with  the  radar  line  of  sight.  As  described  by  Wehner 
in(3),  in  this  case  the  "effective  rotation  axis"  is  coplanar  with  the 
radar  line  of  sight  and  the  true  rotation  axis  and  perpendicular  to  the 
radar  line  of  sight.  The  effective  rotation  vector  has  magnitude  of  the 
vector  projection  of  the  true  rotation  vector.  Also,  the  delay-doppler 
image  plane  is  normal  to  the  effective  rotation  axis. 


2. 


This  chapter  contains  a  dascription  of  two  inpleaen tat ions  of  ISAR 


pr ocas sing .  Thay  diffar  in  tha  waveforms  transmitted.  Tha  first  usas  a 
linaar-FM  pulsa  waveform.  Tha  aacond  uaas  a  stepped- frequency 
pulsa- burst  waveform.  Tha  aacond  waveform  is  a  discretization  of  tha 
linaar-FM  pulsa  wavafora.  It  is  bast  daacribad  in  coaparison  to  tha 
linaar-FM  pulsa,  and  so  both  tachniquas  ara  daacribad  balow  to  allow  a 
batter  overall  view  of  ISAR  tachniquas. 

2.1  LINEAR -FM  ISAR  IMAGING 

Tha  technique  of  using  a  linaar-FM  pulsa  wavafora  for  ISAR  iaaging 
is  described  in  (2).  A  series  of  pulses  is  transaitted  to  produce  an 
image.  Each  pulsa  has  a  rectangular  envelope  and  linear  swept  frequency 
during  tha  pulsa.  Tha  instantaneous  frequency  of  a  linaar-FM  pulsa  is 
illustrated  in  Figure  2.1. 


/I 


Swtpt-f r*que"cy 
Bandwidth 


Instantaneous  Frequency  of  a  Linaar-FM  Pulse 
Figure  2.1 
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The  received  signal  is  processed  in  the  following  way.  It  is  mixed 
with  the  transmitted  signal  and  sampled,  to  produce  a  discrete  time 
signal  which  when  Fourier  transformed  will  yield  the  range  profile  of 
the  target.  An  illustration  of  the  radar  system  components  along  with  a 
comparison  of  the  transmitted  and  received  signals  for  a  point  target 
and  the  mixer  output  is  shown  in  Figure  2.2.  The  mixer  output  signal  is 
sampled  and  the  discrete  Fourier  transform  (OFT)  is  calculated  to  get  a 
sequence  of  complex  numbers.  The  magnitude  values  of  this  complex  number 
sequence  would  be  the  target  range  profile,  but  in  getting  an  ISAR  image 
several  linear-FM  pulses  are  transmitted  and  received,  mixed,  sampled 
and  transformed.  The  complex  range  profiles  resulting  from  each  pulse 
are  aligned  in  rows  of  a  two-dimensional  array  so  the  columns  represent 
a  range  bin,  the  same  range  value.  This  two-dimensional  array  corres¬ 
ponds  to  the  image  plane  shown  in  Figure  1.1  such  that  the  rows,  each 
containing  a  mixer  output  sequence  resulting  from  a  pulse,  are  aligned 
parallel  to  the  y-axis.  The  columns  are  parallel  to  the  x-axis. 
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The  DFT  of  each  column  in  the  array  is  than  calculatad.  This  sorts  the 
doppler  frequencies  of  the  data  in  the  range  bins  yielding  cross-range 
profiles  when  the  magnitudea  of  the  resulting  complex  numbers  are 
calculated.  After  this  second  step,  the  array  cot.  tains  an  image  which 
shows  the  distribution  of  reflected  energy  over  the  surface  of  the 
target . 

2.2  STEPPED -FREQUENCY  ISAR  IMAGING 

This  variation  on  ISAR  techniques,  from  (3),  may  be  described  as 
using  a  discretized  version  of  the  linear-FM  pulse  waveform.  The 
stepped- frequency  waveform  consists  of  a  sequence  of  narrow  pulses,  with 
rectangular  envelopes,  cosq>rising  a  burst  of  pulses.  This  burst  of 
pulses  is  illustrated  in  Figure  2.3. 


•>>•*>•  2.1.  >!■»»»<  Fr»^M»wcy  Bunt 
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Each  succeeding  pulse  in  the  burst  is  stepped  in  frequency.  The 
separation  between  pulses  is  large  enough  for  the  transaitted  pulse  to 
travel  to  the  target  and  the  echoed  pulse  to  return  to  the  receiver 
before  the  next  pulse  is  transaitted.  This  is  illustrated  in  Figure 
2.4,  showing  the  relative  delay  of  an  echo  froa  a  point  target  coopered 
to  the  tiae  interval  between  transaitted  pulses.  The  interpulse  spacing 
is  referred  to  as  tiae  T^.  The  duration  of  the  burst,  referred  to  as 
tiae  T^,  is  the  same  length  as  a  linear-FM  pulse,  yielding  the  same 
doppler  resolution.  This  comparison  is  based  on  the  formula  for 
cross-range  resolution,  froa  (3), 


4sin(0/2)  20 

For  a  linear-FM  pulse  wavefora,  A  is  the  wavelength  of  the  center 
frequency,  whereas  for  a  stepped-frequency  pulse-burst,  A  is  the 
wavelength  of  the  lowest  frequency  is  the  burst.  The  total  angle 
rotated  by  the  target  is  0,  and  target  angular  velocity  is 


14  HataMwi  •«  T 
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In  Figure  2.5,  the  stepped -frequency  burst  of  pulses  is  compared  to  the 
linear-FM  pulse,  showing  the  same  swept -frequency  bandwidth  transmitted 
during  the  respective  waveforms.  This  bandwidth  is  the  same  for  each, 
to  yield  equal  delay  resolution.  The  equation  for  delay  resolution  is 
as  follows,  from  (3), 

where  B,  also  used  in  Figure  2.5,  is  the  swept -frequency  bandwidth  of 
the  linear-FM  pulse  and  the  stepped-frequency  burst  and  c  is  the  speed 
of  light.  It  is  seen  that  the  stepped-frequency  burst  is  a  discretized 
version  of  the  linear-FM  pulse.  Each  pulse  in  the  stepped-frequency 
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waveform  has  duration  ,  chosen  to  correspond  to  the  two-way  travel 
time  of  a  signal  across  the  field  of  view  of  the  image  along  the  delay 
dimension.  This  arrangement  enables  the  following  situation,  as  illu¬ 
strated  in  Figure  2.6.  The  reflected  pulse  arriving  at  the  receiver  can 
be  considered  to  be  composed  of  component  pulses  reflected  by  resolution 
cells  distributed  in  delay.  By  measuring  the  mixed  received  signal  at 
the  instant  of  time  T^/2  seconds  after  the  round-trip  delay  to  the 
center  of  the  image  area,  illustrated  as  x  in  Figure  2.6,  the  mixer 
output,  a  complex  quantity,  will  be  the  sum  of  received  signal  com¬ 
ponents  from  each  resolution  cell  in  the  target.  A  measurement  of  the 
complex  signal  is  made  in  this  fashion  for  each  return  pulse,  after 
mixing  with  the  transmitted  pulse,  resulting  in  a  sequence  of  complex 
numbers  for  each  burst  of  pulses. 


Several  bursts  are  transmitted  to  produce  an  image.  The  sequence  of 
received  measurements  for  each  burst  can  be  placed  in  rows  of  a  two 
dimensional  array  in  the  same  way  as  linear-FM  processing.  The  rows 
represent  the  spectrum  of  a  range  profile  of  the  target  and,  when  the 
inverse  Fourier  transform  is  computed,  result  in  the  complex  range 
profile.  To  generate  an  image,  the  forward  Fourier  transform  magnitude 
of  each  column  is  computed,  as  it  was  for  the  linear-FM  technique,  to 
sort  doppler  frequency  components  into  a  cross-range  profile.  The  result 
is  a  two  dimensional  image. 

There  are  several  advantages  to  using  the  stepped- frequency  tech¬ 
nique  to  make  ISAR  images.  Since  a  narrow  pulse  of  constant  frequency 
is  the  basic  unit  of  the  waveform,  with  each  of  these  separated  by  a 
relatively  long  time  interval,  the  wide  instantaneous  bandwidth  and  high 
sampling  rate  requirements  of  transmit  and  receive  equipment  is  removed 
compared  to  what  is  necessary  for  the  linear-FM  technique.  The 
stepped- frequency  ISAR  technique  is  far  simpler  to  implement  for  this 
reason. 

2.3  ANGULAR  DOPPLER  PROCESSING 

As  described  by  Mensa  (2),  in  the  stepped -frequency  waveform  and 
linear-FM  waveform,  and  associated  imaging  techniques,  there  can  be  seen 
a  duality  of  time  and  angle  of  the  target  during  signal  transmission. 
The  target  rotates  continuously  and  there  is  a  one-to-one  correspondence 
between  time  and  angular  position.  Therefore,  the  Fourier  transforms, 
which  were  before  seen  as  operations  on  temporal  signals,  can  also  be 
viewed  as  transforms  of  spatial  signals. 

Making  use  of  this  concept  means  that  the  stepped-frequency  ISAR 
approach  can  be  used  in  a  laboratory  setting  to  make  two  dimensional 


3.  ALGORITHM  IMPLEMENTATION 


A  computer  program  was  developed  which  implements  stepped-frequency 
ISAR  imaging.  This  was  implemented  as  a  simulation  of  radar  imaging,  so 
targets  were  specified  and  their  scattering  functions  were  determined 
analytically.  The  computer  program  uses  the  scattering  function  to 
calculate  the  mixed  received  signal  from  the  target.  The  ISAR  process¬ 
ing  part  of  the  simulation  uses  the  received  signals  to  generate  an 
image.  This  simulation  was  implemented  in  the  facilities  of  the 
Electronic  Systems  and  Signals  Research  Laboratory  at  Washington 
University.  The  radar  images  obtained  were  displayed  on  a  color 
graphics  display.  To  accomplish  this,  it  was  necessary  to  study  the 
mathematical  model  of  ISAR  processing.  This  model  is  presented  in  this 
chapter.  Also  described  are  computer  algorithms  which  implement  the 
radar  imaging  simulation.  Finally,  the  results  of  the  radar  imaging 
simulation  are  presented  by  examples  of  images  which  were  generated. 

3.1  DISTRIBUTED  TARGET  MODEL 

The  objective  is  a  high  resolution  image  of  the  radar  target.  The 
target  is  distributed  over  many  resolution  cells.  When  the  entire 
target  is  illuminated  by  the  radar,  the  target  is  within  the  antenna 
beam  width,  and  the  back-reflected  signal  is  considered  to  consist  of 
the  superposition  of  the  returns  from  the  multitude  of  resolution  cells. 
Each  resolution  cell  appears  as  a  point  scatterer.  Thus,  the  target  can 
be  viewed  equivalently  as  an  array  of  point  scatterers  which  produces 
the  same  radar  return  as  the  actual  target.  For  the  disk  and  spherical 
targets,  the  resolution  cells  are  located  on  the  surface  as  shown  in 


Figure  3.1  where  one  such  resolution  cell  is  located  as  an  intersection 
in  delay  and  doppler.  In  the  front  view  of  the  target  shown  in  Figure 
3.1,  the  radar  illuminates  the  target  from  a  position  perpendicular  to 
the  plane  of  the  paper. 


Disk  Sphere 


Figure  3.1  The  Disk  and  Sphere  Targets 

In  fact,  for  the  spherical  target,  there  are  two  area  elements  on  the 
surface  with  the  same  delay-doppler  coordinates.  This  is  an  ambiguity 
resulting  from  the  geometry  of  the  situation  and  the  properties  of  high 
resolution  radar.  The  radar  resolves  the  three  dimensional  target 
surface  in  only  two  dimensions,  delay  and  doppler.  The  total  reflection 
of  the  resolution  cell  is  the  sum  of  the  reflection  of  the  two  area 
elements. 
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Figure  3  .2.  The  Oeiey-Ooppler  Plane 


The  two  area  elements  of  the  resolution  cell  are  equivalently  mapped 
onto  the  delay-doppler  plane  shown  on  Figure  3.2.  Point  scatterers  on 
the  delay-doppler  plane  can  represent  the  resolution  cells  mapped  onto 

i 

it.  The  reflectivity  of  a  point  equals  the  vector  sum  of  reflectivity 

within  the  resolution  cell,  as  described  in  (3). 

The  two-dimensional  planar  array  of  points  of  reflectivity  squared 

magnitude  is  called  the  scattering  function  of  the  target.  The  signal 

which  would  be  received  from  an  illumination  of  the  point  scatterers  in 

the  scattering  function  is  equivalent  to  the  back  scattered  signal  from 

the  actual  target.  However,  the  following  assumptions  are  imposed  on  the 

reflectivity  of  the  points  during  ISAR  imaging,  as  described  in  (A).  It 

is  assumed  the  reflectivity  does  not  change  over  the  bandwidth  of 

illuminating  signals.  Also,  it  is  assumed  the  reflectivity  does  not 

change  as  target  aspect  changes  with  rotation  during  the  sequence  of 
I 
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3.2  STEPPED- FREQUENCY  ISAS  PROCESSING  MODEL 

At  the  instant  of  mixer  output  measurement  specified  as  T.  in 

b 

Figure  2.6,  the  received  signal  is  comprised  of  the  superposition  of 
components  echoed  from  the  point  scatterers  in  the  target.  The  super¬ 
position  occurs  over  the  delay-doppler  plane  as  follows.  The  whole 
received  signal, 

/v  .  .  a;  .  . 

s  (t)  -Is  (t), 

1  u,v  r 

where  u  is  the  delay  location,  v  is  the  doppler  location,  and  t  is  the 
time. 

The  received  signal  component  due  to  a  resolution  cell  reflection  is, 

sr(t)  -  b  sT(t-r), 

Aj 

where  t  is  the  travel  time  delay,  and  b  is  reflectivity  at  the 
delay-doppler  coordinate.  This  received  signal  is  a  reflection  of 
the  transmitted  signal, 

~  , .  x  'Z , .  ,  j  2rf  .  t 

sT(t)  ■  f(t)  eJ  l  , 

VneTe  i  is  the  pulse  number,  fj  is  the  frequency  of  the  pulse,  and 
f(t)  is  the  complex  envelope  of  the  transmitted  pulse.  Assume  f(t)  * 

1.0  within  the  pulse. 

For  the  discussion  in  this  section,  the  reflectivity  is  assumed  to 
be  a  real  number  and  a  deterministic  value  equal  to  the  square  root  of 
>(u,v),  reflected  power  in  the  resolution  cell  at  the  delay-doppler 
coordinate  (u,v).  The  quantity  A(u,v)  is  also  called  the  scattering 
function.  In  this  section,  the  processing  is  shown  to  map  correctly  the 
delay-doppler  reflected  power,  a  non-randomly  fluctuating  quantity. 
This  will  verify  the  proper  operation  of  the  processing. 


If  ISAS  processing  works  correctly,  it  reconstructs  froa  the 
received  signals  the  scattering  function  which  produced  then,  Because 
of  the  linear  superposition  property  of  the  Fourier  transfora,  the 
iaaging  of  the  entire  scattering  function  is  the  sta  of  the  iaeging  of 
single  point  scatterers.  Therefore,  consider  a  single  point  scatterer 
and  the  received  signal  froa  a  pulae  reflected  froa  it, 

SjU)  -  b  •i  b  is  real. 

This  signal  is  a  complex  quantity  consisting  of  real  and  iaaginary, 
alternatively  in-phase  and  quadrature,  coaponents.  To  get  an  1SAR  image 
of  the  single  point,  the  first  step  is  the  mixing  of  the  received  signal 
with  the  transmitted  signal.  This  is  accomplished  by  the  system  shown 
in  Figure  3.3. 
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The  aixer  output  is  derived  «s  follows.  At  th«  instant  of 
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ment  tha  aixer  output  signal 


» 


n>j(t)  *  2*LP(Re(s  j(  t ) )  •Re(s-j'(  t ))  j  <-  j  2*LP[Im(s,i(t))-Im(q(t))] , 


where  LP  means  a  low-pass  filter,  si  is  the  received  signal, 
St  is  the  transmitted  signal,  and  q  ■  ST  shifted  90*  *  j  St • 


In  the  first  tera, 

Re[s.  (t)3  •Re[s_(t)]  *  b  cos  2*f  .  (t-xVcos  2irf.t 

i  T  l  i 

where  Ref  ]  means  the  real  part  is  taken.  Since 

cos(at)*cos(bt)  »  ^[cos(a+b)t  *  cos  (a-b)t], 

it  follows  that 

Re  [  s  t )  ]  *  Re  { Mj( t ) ]  •  ^[cos  2wfi(2t-x)  +  cos  2x^x1. 

Similarly,  since  sin  (at)  .cos  (bt)  ■  ^  [sin(a+b)t  +  sin(a-b)t), 
in  the  imaginary  tera  of  the  mixer  output  signal, 

Ials^t)]  'Ia[q^(t>]  ■  b  sin  2irf  ^t-x)  •  ( -cos  2-irf^t), 

■  ^  l-sin  2irf^(2t-x)  -  sin  2irf^(-x)], 

■  ^  {-sin  2irfi(2t-T)  ♦  sin  2irfix], 


where  Ia[  ]  means  the  iaag inary  part  is  taken.  After  low-pass 
filtering,  the  real  component  of  the  mixer  output  signal  is  given  by 
2*LP[Re(s^(t)) 'Re(q^(t))]  ■  b  cos  2irf^x. 

Also,  after  low-pass  filtering,  the  imaginary  component  of  the  mixer 
output  signal  is  given  by 

j  2*LP[Ia(si(t))>Im(qi(t))]  •  j  b  sin  2*^1. 


Therefor#,  the  total  mixer  output  signal  is  given  by 


m^(t)  *  b  cos  2irf^r  +  j  b  sin  2irf^t, 

or, 

m^t)  ■  b  eJ  i  . 

Tor  a  single  point  target  of  arbitrary  location  in  the  delay- 

doppler  field  of  view,  the  round-trip  delay  t  is  given  by 

t  »  t (t )  -  ”(r-v  t) , 
c  r 

where  vr  is  relative  velocity  of  the  point  target,  and  r  is  initial 
range  to  target. 

There  are  a  fixed  number  of  pulses  in  a  burst,  yielding  a  sequence 
of  mixer  output  values.  Time  within  the  burst  is  incremented  as 
follows: 

t  -  i  T2  +  Tj/2  +  k  T3. 

where  k  is  Burst  index,  Ti  is  pulse  width,  and  i  is  the  pulse  index. 

Using  this  equation  for  time  in  the  above  equation  for  round-trip  delay, 
t  may  be  expressed  as 

t  -  |(r  -  vr(i  T2  ♦  Tj/2  ♦  k  T3)). 

The  expression  for  the  mixer  output  signal  then  becomes, 

■4  -  b  exp[ j  2sfi  J(r  -  vr(i  T2  +  Tj/2  ♦  k  T3>)], 
where  exp  [  ]  means  an  exponential  function. 
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Th#  sequence  of  mixer  outputs  from  e  burst  is  inverse  Fourier 
trensformed  to  get  a  range  profile  RP(n)  of  the  target. 


.  N-l  2x . 

RP(n)  -  IDFTfm. ]  -  5  I  m.  eJ  N  n* 
1  N  i-0  1 


where  N  is  the  number  of  pulses,  n  is  the  range  profile  index,  IDFT 
means  an  inverse  discrete  Fourier  transform  is  coaiputed.  Substituting 

the  expression  for  the  mixer  output  signal  in  this  range  profile 
equation,  it  is  seen  that 

N-i  ,  , 

RP(n)  ■  -  I  exp[j  2wfi  ~  (r-v(i  T2  ♦  T^/2  +  k  T3))]exp[j  —in], 


or  after  rearranging  terms. 


RP(n)  -  £  I  explj  ^  (i  n  +  Hf.  |(r-v  (i  T2  +  Tj/2  +  k  T3>)>} 
i-0 


The  frequency  of  any  pulse  is  expressed  as,  fi  -  fo  +  i  Af,  where  fo 
equals  initial  frequency  used  in  the  burst  and  Af  is  the  frequency  step 

of  the  succeeding  pulses.  It  follows  that 

N_  2 

RP(n)  -  £  I  exp[ j  ^  (i  n  +  — (f  +  iAf)  (r  -  V  (i  T-  +  T./2  +  kT,))) 
N  N  CO  r  L  1  J 

2B  1 

By  using  NAf  -  B,  and  —  -  — ,  the  above  expression  for  the  range 
profile  is  simplified.  Also,  the  following  restriction  on  point  target 


motion  is  imposed, 


|vr(iT2+T1/2  +  kT3)|<Ar, 
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where  Ar  is  rang*  resolution.  This  restriction  implies  that  the  point 
does  not  move  out  of  its  resolution  cell  during  rotation.  This  results 
in  the  following  expression  for  the  range  profile, 

N-l 


RP(n)  -  £  I  exp[j  jj”  (i  n  ♦  ~  -f  Nf0  “)]. 


i-0 


N 


Ar 


(3.1) 


Since  —  “ 
Ar 


-  n  ,  and  —  -  t.  ,  it  follows  that 

O  C  K 


.  b  j  2*f  t  j  —  i(n-n  ) 

RP(n)  ■  r  eJ  ok  l  eJ  N  o 

N  i-0 


where  r  is  a  delay  time  associated  with  a  particular  burst.  The 
k 

resulting  range  profile  is  then  froa  (3), 


(3.2) 


.  .  i  2trf  t  sin  *(n-n0) 

RP(n)  -  b  eJ  Z  o  k  - — 


sin  — (n-nQ) 

The  range  profile  sequence  has  a  function  narrow  "impulse- like" 
term  which  has  peak  response  within  the  down-range  resolution  cell  where 
the  point  target  is  located,  low-level  sidelobes  next  to  the  peak  and 
essentially  zero  response  in  other  positions.  Each  burst  has  a  range 
profile  described  by  the  above  equation,  containing  a  phase  term  which 
is  a  function  of  burst  number  and  target  velocity.  In  the  following, 
the  above  equation  is  simplified  by  replacing  the  narrow  "impulse- like" 
function  term  with  an  ideal  impulse  which  is  unity  within  the  down-range 
resolution  cell  of  the  point  target  and  zero  in  other  positions,  thus 
ignoring  the  low-level  sidelobes  of  the  above  response. 

In  this  manner,  the  above  equation  is  expressed  as, 

RP(n)  -  b  ej  2vf°  Tk  6(n-n0)  (3.3) 


m 
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In  the  next  step  of  ISAR  processing,  the  point  target  is  located  in 
cross-range  position  by  taking  the  discrete  Fourier  transform  (DFT)  of 
the  sequence  of  range  profile  values  of  constant  range, 

DFT  [RP(no>]  -lb  exp[.j  2xfoTkJ  exp[-j  jp  km}. 

''  k-0 

This  is  described  as  a  cross-range  profile, 

N-l  - 

CRP(m)  -  b  I  exp[-j  rr  ^mk  "  Nf  t  )], 

I  a  N  OK 

k»0 

2 

where  m  is  the  doppler  index  of  the  sequence  and  *  -  (r-V  k  T^). 

Using  this  expression  for  burst  delay  time  in  the  cross-range  profile 
results  in 


CRP(m)  -  b  I  exp! -j  {—  (km  -  Nf0  —  +  kT3  J  VrNf0)]. 

k-0 


fc  -  — ,  the  cross-range  profile  becomes. 


CRP(m)  -  b  I  exp[-j  ~  k(m  +  NT.,  7  V)]  exp[j  2xf  t ] . 

k-0  N  3  A  r  o 


N*1  2  NT-2 

CRP(m)  -  b  exp[j  2irf  t]  Z  exp[-j  Tp  k(m  -  — —  m  Ax  u>  )]. 

0  k-0  N  A  o  r 


i  r. *  k( m-m  )  .  1  2ef  t  ,sin  "  (n_no), 

J  N  o  -  be  [ - J. 


sin  -  (n-nD) 
N 
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By  retaining  only  the  peak  impulse  response  and  ignoring  low-level 
sidelobes, 

CRP(n)  ■  b  e^  ^^o  6(m-m  ). 

o 

The  magnitude  of  the  cross-range  profiles  forms  the  image.  When  simpli- 
f lying  the  expression  by  ignoring  side- lobes  and  using  only  the  Ideal 
impulse  peak  response,  the  profile  magnitudes  form  an  image, 

|CRP(m)|  -  Image  (n,m)  ■  b  5(n-n0,  m-mo). 

This  is  the  final  result  of  ISAR  processing.  It  is  seen  that  the 
point  target  shows  up  in  the  image  at  delay-doppler  coordinate  (nQ,  nO 
as  a  point  with  magnitude  equal  the  square  root  of  the  scattering 
function,  as  it  is  supposed  to. 

For  some  parameter  choices,  the  total  angle  of  target  rotation 
during  the  image  frame  time  is  more  than  a  few  degrees.  In  this  case, 
the  point  target  described  above  will  migrate  out  of  its  resolution  cell 
during  the  image  frame  time,  invalidating  the  assumption  used  in  getting 
equation  3.1  above.  The  effect  will  be  a  distorted  image.  This  problem 
is  described  in  (3)  and  is  known  as  "range  walk".  A  resulting  image 
will  appear  as  if  each  range  profile  is  staggered  in  position  relative 
to  its  neighbor.  This  distortion  is  corrected  by  multiplying  the  mixer 
output  samples  by  a  correction  factor  derived  here. 

Compensating  for  the  range  walk  problem  results  in  the  following 
modification  of  equation  3.1, 

b  N-i  2  it  i  2NfQ 

RP(n)  •  -  I  exp[ j  —  (i  n  +  (r  -  VkT_)  +  - -  (r  -  VkT_))]. 

N  N  AT  J  C  J 
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This  range  profile  equation  can  be  reduced  to, 


RP(n)  -  jj  exptj  ^  fQTk  ]  ^  exp[ j  —  i((n-nQ)  -  k  ^)].  (3. A) 

By  representing  the  response  as  an  ideal  impulse  and  ignoring  low-level 

sidelobes,  this  range  profile  is  expressed  by  the  following  function 

additionally  dependent  on  burst  number, 

b  VT3 

RP(n,  k)  -  g  6  (n-no  -  k  -^). 

This  expression  for  the  range  profile  looks  like  that  in  equation 
3.3,  on  page  22,  except  that  for  each  succeeding  burst,  the  profile 
position  is  offset  in  range  position  by  a  term  due  to  target  motion 
between  bursts.  This  is  corrected  by  multiplying  the  mixer  output 
sequence  for  a  burst  by  a  correcting  factor  prior  to  computing  the 
transform.  This  modifies  equation  3.4  so  that  it  looks  like  equation 
3.2,  resulting  in  a  corrected  range  profile, 

RP(n)  -  jj  exp[j  ^  Tfc]  V  exp[j  ^  i(n-no  -  k  ^)]S.k> 

where  is  the  range  walk  correction  which  equals 

2  it  VT3 

exptj  jj-  ik  ~  ].  Then, 

b  2ir  N- 1  2tt  VT3  VT3 

RP(n)  -  jj  exptj  jj-  f0xk]  ,£o  exptj  ~  i((n-n  )  -  k  ~)]exp[j  J-  ik 

b  ,  .  2it  .  N-l  .  2ir  . ,  . . 

“  jj  exptj  jj-  f0TkJ  exptj  jj-  i(n-n0)J. 

From  this,  it  is  seen  that  the  range  walk  correction  factor  does  restore 
the  range  profile  to  the  proper  form. 


3.3  COMPUTER  SIMULATION 

A  FORTRAN  program  was  written  to  achieve  the  goal  of  producing  ISAR 
images  in  the  laboratory  by  simulation.  The  dimensions  of  simulated 
targets  were  specified,  and  their  scattering  functions  were  supplied  for 
use.  The  simulation  contained  two  parts.  The  first  part  inputs  the 
scattering  function  and  calculates  the  mixer  output  signals  from  that 
target  for  all  pulses  and  bursts  transmitted.  The  second  part  processes 
the  signals  using  Fourier  transform  techniques,  as  described  here  ,  to 
produce  an  image  data  file.  An  image  display  routine  is  used  to  display 
the  simulated  ISAR  image  on  a  color  monitor. 

An  algorithm  is  presented  in  Figure  3.4  to  describe  the  programming 
of  the  part  of  the  simulation  which  calculates  mixer  output  signals. 


Algorithm  of  Mixer  Output  Simulation 
Figure  3.4 


A  two-dimensional  array  is  used  to  store  the  resulting  radar  signals; 
each  burst  forms  a  row  of  the  array. 

Two  analytic  targets  were  chosen,  a  rough  disk  and  a  rough  sphere. 
The  field-of-view  of  the  image  was  chosen  to  be  76.8  meters  in  both 
range  and  cross-range  dimensions.  The  resolution  in  both  dimensions  was 
chosen  to  be  60  centimeters.  The  scattering  function  is  a  128  by  128 
array  of  reflected  power  values.  The  diameters  of  both  the  disk  and 
sphere  were  chosen  to  be  38.4  meters  with  both  centered  in  the  field- 
of-view.  It  was  desired  to  write  a  computer  program  that  would  generate 
an  image  of  size  256  pixels  on  each  side  with  the  same  resolution  of  60 
cm  in  each  coordinate.  This  is  accommodated  in  the  Fourier  transform 
techniques  by  padding  with  zeroes.  This  results  in  an  image  that  is 
magnified  by  a  factor  of  two  and  is  easier  to  see  on  the  color  monitor. 

In  this  simulation,  128  pulses  per  burst  and  128  bursts  are  used  in 
the  transmitted  waveform.  This  covers  the  field-of-view  of  interest, 
given  the  resolution.  The  resulting  two-dimensional  array  holding  the 
mixer  output  signal  is  of  size  128  by  128.  The  range  profiles  are 
calculated  with  an  inverse  discrete  Fourier  transform  (DFT).  This 
procedure  uses  a  subroutine  for  DFT  calculation  which  was  supplied  by 
the  University.  Also  in  computing  a  range  profile,  the  128  point  signal 
is  placed  in  the  center  of  a  256  point  array,  padded  on  both  ends  with 
64  zeroes,  for  the  inverse  DFT  operation.  The  resulting  256  point  range 
profile,  divided  by  the  number  of  points,  replaces  the  row  of  signal 
data  in  the  two-dimensional  array. 
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In  the  final  step  of  calculations,  the  cross-range  profiles  of  the 
image  are  made  using  a  forward  DFT  on  the  columns  of  the  two-dimensional 
array  resulting  from  the  last  operation.  However,  rearrangement  of  the 
data  in  the  256  point  array  is  necessary.  Instead  of  placing  the  128 
data  points  in  the  center  of  the  256  point  array,  the  first  64  data 
points  are  put  in  the  first  64  array  locations  and  the  last  64  data 
points  are  placed  in  the  last  64  array  locations,  leaving  the  center 
array  locations  zeroed.  This  is  due  to  the  way  the  DFT  routine 
operates.  After  the  DFT  is  computed,  the  zero  frequency  location  is  in 
the  first  array  location,  so  the  data  are  rearranged  with  the  last  half 
of  the  256  points  switched  with  the  first  half.  This  is  illustrated  in 
Figure  3.5.  The  squared  magnitudes  of  the  resulting  complex  numbers  are 
stored  in  each  column  of  the  two-dimensional  array.  The  resulting 
two-dimensional  image  data  array  is  used  with  a  color  scale  to  display 
the  image . 
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The  computer  listing  is  included  in  the  Appendix  in  Figure  7.1.  In 
the  main  routine,  PR0C3,  the  target's  scattering  function  is  input  and 
the  mixer  output  signal  is  calculated.  This  signal  is  normalized  to 
peak  magnitude.  Subroutine  IMAGE2D  is  then  called  to  reconstruct  the 

image. 

As  a  part  of  this  simulation  development,  a  modification  of  the 
Fourier  transform  techniques  described  above  was  investigated.  A 
modified  IMAGE2D  subroutine  was  constructed  in  which  the  Fourier  trans¬ 
forms  were  "windowed"  with  a  normal  function, 

w(n)  «  e  -  0.09  (n-N/2)2  #  n  „  index,  N  ■  128. 

Windowing  is  described  as  multiplying  the  complex  data  sequence  by 
the  window  function  before  calculating  the  transform.  Due  to  the 
synsetry  properties  of  the  Fourier  transform,  each  of  the  delay  and 
doppler  transforms  are  windowed  with  the  same  function.  This  modified 
routine  is  called  IMAGE2DW  and  is  listed  in  the  appendix. 
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COMPUTER  FACILITIES 

The  computer  program  was  written  in  FORTRAN  or  a  Massachuse* t s 
Computer  Company  (MASSCOMF'  image  processing  computer,  modei  MC5  3-'.  .  with 
a  color  graphics  display  terminal.  The  color  of  the  pixels  m  the 
raster  displav  are  controlled  by  a  color  scale  with  64  intensity  levels. 
A  heated  object  color  scale  was  chosen.  Once  a  two-dimensional  arr,v  of 
data  is  available,  an  image  display  routine  is  executed  which  reads  the 
data  and  generates  the  image  or  the  screen.  A  picture  of  this  facil;*.' 
is  presented  in  Figure  3.6. 


□ 


Fhetograph  of  Computer  facilities 


Figure  1.6 
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3.5  SIMULATION  RESULTS 

Several  scattering  functions  with  the  same  resolution  and  size  were 
used  in  tests  of  the  ISAR  imaging  simulation.  These  scattering  func¬ 
tions  can  be  displayed  to  see  what  the  actual  target  looks  like.  Their 
data  arrays  are  sized  128  by  128,  but  it  is  desired  to  see  them  as  a  256 
by  256  pixel  image.  Therefore  one  data  point  drives  4  pixels  in  order 
to  double  the  size  of  the  image.  As  first  test  cases,  point  target 
pairs  were  used.  One  pair  has  a  point  in  the  center  of  the  image,  at 
location  (64,64)  of  the  128  by  128  scattering  function,  and  a  point  at 
(64,96).  Another  point  target  pair  uses  a  center  point  and  a  point  at 
(64,66).  The  ISAR  imaging  simulation  was  used  to  reconstruct  ISAR 
images  of  the  scattering  functions.  These  radar  images  are  displayed  in 
Figure  3.7  with  the  closely  spaced  pair  on  the  left  (a).  The  resulting 
images  are  256  by  256  pixels  in  size  due  to  the  processing  implemen¬ 
tation  described  above.  The  fact  that  the  radar  images  look  like  the 
scattering  functions  confirms  the  proper  operation  of  the  ISAR  process¬ 
ing  for  these  test  targets. 


(a)  (b) 

ISAR  Images  of  Point  Target  Pairs 
Figure  3.7 
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During  image  display,  the  image  data  is  scaled  to  the  color  scale 
so  that  the  peak  value  is  shown  as  white  and  zero  values  appear  as 
black.  In  the  images  in  Figure  3.7,  the  peaks  of  the  responses  from  the 
point  targets  are  seen  very  well.  However,  there  are  also  sidelobes  to 
the  peaks  which  are  subdued  in  these  pictures  by  the  color  scale  and  the 
effects  of  film  processing. 

The  scattering  functions  of  the  rough  disk  and  rough  sphere  are 
displayed  in  Figure  3.8.  These  scattering  functions  come  from  reference 
(1).  The  disk  is  shown  on  the  left  (a)  in  the  figure.  Note  that  the 
sphere  (b)  shows  the  effects  of  shadowing  which  makes  the  illumination 
look  somewhat  like  a  partial  moon.  Also,  the  reflected  power  is 
greatest  from  the  point  on  the  target  closest  to  the  radar  and  decreases 
for  points  closer  to  the  shadowed  edge.  In  this  figure  and  in  other 
similar  images  shown  here,  the  radar  is  illuminating  from  the  left  side. 


Scattering  Functions  of  Disk  and  Sphere 
Figure  3.8 
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Reconstructed  images  of  the  disk  und  sphere  are  presented  in  Figure 
3.9.  Horizontal  corresponds  to  the  delay  dimension  and  vertical  to  the 


doppler  dimension.  Note  that  the  sphere's  image  (b)  looks  like  that  in 
Figure  3.2.  The  disk's  image  is  on  the  left  (a)  in  Figure  3.9.  There 
is  substantial  low-level  sidelobe  energy  in  the  image  which  can  be 
smoothed  out  by  windowing,  but  results  in  a  loss  of  resolution. 


(a)  (b) 

ISAS  Images  of  Disk  and  Sphere 
Figure  3.9 


An  image  of  a  center  located  point  target,  reconstructed  using 
windowed  transforms  is  shown  in  Figure  3.10.  The  window  is  described  at 
the  end  of  Section  3.3.  As  can  be  seen,  windowing  blurs  a  point  target. 
The  profile  of  the  image  is  also  plotted  in  the  figure.  The  plot  also 
shows  the  form  of  the  window.  A  similarly  windowed  image  of  the  spheri¬ 
cal  target  is  presented  in  Figure  3.11  along  with  the  profile  plot. 
Those  images  suggest  that  the  blurring  introduced  by  windowing  reduces 
image  quality,  but  side-lobe  energy  is  reduced. 
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Windowed  ISAR  Image  of  Spherical  Target 
and  Profile  Plot 


Figure  3.11 


4.  STATISTICAL  TARGET  MODEL 


In  the  preceding  chapter,  the  reflectivity  of  the  target  surface, 
for  each  resolution  cell,  was  taken  to  be  a  fixed  real  value  equal  to 
the  square  root  of  the  scattering  function  value.  This  was  for  the 
purpose  of  verifying  the  ISAR  processing.  The  next  desired  step  was  to 
image  a  diffuse  target  surface.  The  reflectivity  of  a  diffuse  target  is 
a  complex  random  process. 

The  resolution  cells  are  distributed  in  the  delay-doppler  plane. 
Each  is  described  as  a  slowly  fluctuating  point  target  in  the  presence 
of  additive  noise.  The  resolution  cell  is  made  up  of  several  reflecting 
surfaces  which  combine  to  produce  a  random  reflected  signal.  Slowly 
fluctuating  means  the  reflectivity  is  partially  correlated  relative  to 
the  transmitted  signal  pulse  burst  period. 

4.1  MODEL  OF  A  SLOWLY  FLUCTUATING  POINT  TARGET 

The  source  for  this  model  is  Van  Trees  (5).  The  radar's  trans¬ 
mitted  signal  is  represented  as  follows: 

sT(t)  »  A  f(t)  exptjZir^t] 

mi 

where  f(t)  is  the  normalized  complex  envelope  of  transmitted  signal, 

A  is  the  amplitude,  and  f^  is  the  carrier  frequency.  The  power  of  such 
2 

a  signal  is  A  ,  from  (6,7). 


Consider  a  point  target  in  space,  representing  a  resolution  cell  in 
the  delay-doppler  plane  with  zero  velocity.  The  reflected  signal  from 
such  a  target  includes  the  effect  of  the  random  reflectivity  process. 
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This  noise  source  is  Multiplicative  and  the  received  signal  can  repre¬ 
sented  as  follows: 

sr(t)  -  b  f (t-r )  exp[j  2*f ^(t-T ) ] , 

A* 

where  t  is  round  trip  delay.  The  tern  b  is  a  complex  Gaussian  random 
variable  with  sero  mean.  The  variance  of  the  real  and  imaginary  com¬ 
ponents  is  equal  to  half  the  power  of  the  reflected  signal.  A  point 
target  with  non-sero  velocity  will  have  the  same  form  of  reflected 
signal,  but  the  carrier  frequency  will  be  shifted  by  the  Doppler  effect. 

A  target  which  is  extended  in  the  delay-doppler  plane,  and  slowly 
fluctuating,  is  considered  to  be  composed  of  an  array  of  point  targets 
representing  the  resolution  cells,  each  contributing  a  reflected  signal 
of  the  above  form.  So,  the  total  received  signal  from  such  a  target  is 
the  sum  of  the  received  sign^b  from  each  point  target.  Representing 
this  sum  and  the  additional  fecture  that  the  signal  has  a  doppler 
frequency  shift,  the  received  signal  can  be  expressed  as, 

sr(t)  -Zb  f(t-Tu)  exp[j(u)i  +  u.r)t], 
u,v 


where 


XV 

_ r  is  the  doppler  frequency,  -  2*^,  the  transmitted 


frequency,  and  (u,v)  is  an  element  of  the  delay-doppler  plane. 


The  random  reflectivity  term  has  the  property  that  it  is  wide-sense 
stationary  and  spatially  uncorrelated.  This  property  along  with  the 
slowly  fluctuating  property  implies  an  approach  in  simulating  the  mixer 
output  from  a  target. 

For  a  resolution  cell,  the  reflectivity  from  pulse  to  pulse  is 
correlated  by  using  a  recursive  filter.  White  Gaussian  noise  is 
filtered  so  it  has  a  desired  autocorrelation  function.  Two  filters  are 
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operated  to  generate  separate  outputs  for  aach  of  tha  raal  and  imaginary 
comp^.._nt*  of  tha  raf lactivity. 

Tha  following  racursiva  filtar  was  chosan  from  (8), 

Y„  *  «Y„  ,  ♦  (I"*2)172  *  • 

n  n- 1  n 

Tha  input  par asm  tar  X  is  a  whita  gauasian  variable  which  has  a  waightad 

n 

contribution  to  tha  output  valua  Y^.  Tha  pravious  Y^_  ^  valua  has  a 
waightad  contribution  to  tha  output.  Tha  autocorralation  function  of 
tha  filtar  output, 

Ry(k)  »  Rx(o)  alkf  ■  ox2  al^l  . 

Tha  shapa  of  tha  autocorralation  function  is  shown  in  Figure  4.1.  The 
slopa  of  tha  function  naar  tha  peak  is  used  to  ralata  the  parameter  a  to 
correlation  interval  parameter  Cl  . 
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ISAR  processing  calculates  the  Fourier  transform  of  the  delay- 
doppler  signal  having  random  reflectivity.  The  result  can  be  viewed  as 
a  convolution  of  the  desired  power  density  of  the  target  with  a  poor 
estimate  of  the  power  spectrum  of  the  random  variable,  as  described  in 
(6,  Chapter  11,  Section  3).  This  signal  is  itself  a  random  process. 
The  power  spectrum  of  the  correlated  reflectivity  components  is  given  in 
(7,8,9), 

Sy(k)  -  |H(k) | 2  Sx(k), 

where  S^Ck)  is  the  power  spectrum  of  the  white  Gaussian  noise,  Sy(k)  is 
the  power  spectrum  of  the  filtered  noise  and  H(k)  is  the  filter's 
transfer  function.  However,  the  Fourier  transform  of  a  limited  sample 
of  a  random  process  is  itself  a  randomly  varying  quantity.  So,  the 
convolution  results  in  random  signal  and  resulting  ISAR  images  are 
expected  to  look  noisy. 

The  effects  of  additive  noise  are  also  to  be  included  in  the 
simulation  of  the  received  signal.  Adding  this  term  to  the  equation  for 
the  received  signal  yields  the  following: 

S  (t)  -  I  f(t-t)  bej  2’fi(t'T)  +S(t) 

u,v 

The  last  term  is  a  complex  gauss ian  random  process  with  zero  mean. 
It  is  included  to  model  wideband  channel  and  receiver  noise.  The  first 
part  of  the  signal  is  normalized  to  power  equal  1.0,  and  the  variance  of 
the  additive  noise  term  is  varied  to  produce  a  desired  signal-to-noise 
power  ratio  (SNR), 

P 

SNR  -  —  . 
n 
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where  P  and  P  are  signal  and  noise  power  over  the  receiver  bandwidth, 
s  n 

P  A2 
s  A 

Further,  SNR  *  p—  *  — r  ,  where  A  is  the  signal  amplitude,  and  since 
n 

A2  -  1,  then  B  *  1/ JSNR.  Thus  the  standard  deviation  of  the  real  and 
imaginary  terms  of  complex  white  Gaussian  noise  are  B//7. 

Because  of  the  linearity  of  the  Fourier  transform,  the  result  of 
ISAR  processing  the  above  mixer  output,  with  the  additive  noise,  will 
add  to  the  image  a  signal  which  is  the  Fourier  transform  of  the  sampled 
white  gaussian  noise.  Since  the  noise  is  white,  it  is  expected  the 
image  will  have  noise  added  to  it  all  over  the  field  of  view. 


ft. 


4.2  IMPLEMENTATION  IN  THE  SIMULATION 


The  computer  simulation,  written  in  FORTRAN,  was  modified  to 
calculate  a  received  mixer  output  signal  from  a  diffuse  target.  The 
scattering  function  of  the  target  is  input  into  the  simulation,  which  is 
a  delay-doppler  mapping  of  the  reflected  signal  power  for  each 
resolution  cell.  The  target  is  slowly  fluctuating,  so  the  random 
reflectivity  of  each  call  of  the  scattering  function  is  statistically 


correlated  in  time.  This  implies  an  outer  loop  which  increments 
resolution  cell.  The  algorithm  of  this  computer  program  is  presented  in 
Figure  4.2. 


a 


Read  in  Scattering  Function 

Find  Sum  of  Total  Power  within  Scattering  Function 
Initialize  20  Array  for  Mixer  Outputs  to  Zero  (FRAME) 

Increment  Resolution  Cell  ( i .  J) 

Initialize  Random  Reflectivity  Recursive  Filter 
Increment  Burst  Number  m 
Increment  Pulse  Number  n 

Pick  Correlated  Random  Reflectivity 
Calculate  Mixer  Output  S(n) 

Sum  Mixer  Output  to  2D  Frame  (n.m) 

End  Loop 
End  Loop 
End  Loop 

Calculate  Desired  Noise  Standard  Deviation  from  SNR 
Increment  Burst  Number  n 
Increment  Pulse  Number  n 

Pick  White  Gaussian  Noise 
Normalize  Mixer  Output 
Add  Noise  to  Mixer  Output 
End  Loop 
End  Loop 
End 


Algorithm  of  Mixer  Output  Simulation 
Using  a  Correlated  Statistical  Target 
and  Additive  White  Noise 
Figure  4.2 
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For  each  resolution  cell  in  the  scattering  function,  one  at  a  time, 
the  mixer  output  for  all  pulses  in  all  bursts  is  calculated  using  random 
number  generators  to  get  a  complex  Gaussian  random  number  for  reflec¬ 
tivity.  The  mixer  outputs  are  stored  in  a  two- dimensional  array.  The 
whole  received  signal  is  the  point-to-point  sum  of  two-dimensional 
arrays  from  each  cell. 

After  all  cells  have  been  incremented,  and  the  two  dimensional 
array  contains  the  mixer  output  resulting  from  the  whole  target,  the 
mixer  outputs  are  normalized  to  power  equal  1.0.  Then  Gaussian  random 
number  generators  are  used  for  additive  noise  with  power  adjusted  to  get 
a  desired  signal-to-noise  ratio. 

Returning  to  the  part  of  the  algorithm  where  the  mixer  output  for  a 
pulse  is  calculated,  the  mathematics  which  is  used  to  carry  this  out  is 
developed  in  the  following  way.  As  described  before,  the  received 
signal  is  composed  of  components  from  all  resolution  cells  in  the  target 

A/ 

and  includes  added  noise.  Assuming  the  received  pulse's  envelope  f(t-x) 
is  unity,  the  received  signal  can  be  expressed  as, 

s  (t)  ■  I  b  eJV  i  r  +  n(t). 
r  u,v 


This  signal  s  (t)  is  mixed  with  e^  ^i  t,  the  carrier  signal,  to  give  the 


mixer  output  signal  which  can  be  written  as. 


s v  ^  1  Ui  t 

s  ■  Z  b  e  J  r  +n. 


m 


u,v 


m 


Letting  9 ,  ■  wft,  with  i  indicating  pulse  'number ,  and 


b  •  b(u,v)  ■  |b|e^b  then, 


m 


Z  fb|ej(  VV  +  n 


u,v 


Since  cos($.-8.)  ■  cos$,  cos0.  +  sin$.  sin0.  , 
o  .  i  o  1  D  1 


and  sin($,-0.)  *  sin$.  cos0.  -  cos$.  sin0.  , 
b  i  b  i  b  l 


s  ■  I  {(|b|  cos$.  cos©.  +  }b|sin<».  sin8.) 
m  b  i  b  l 


U,v 


+  j(|b^sin$,  ccs0 .  -  |b|cos<>,  sin6.)} 
b  l  b  i 


Let  A  *  )b|cos$^,  and  B  *  |b{sin$^,  where  A  and  B  are  Gaussian  random 

2  2  1 

variables  with  zero  mean  and  variance  o,  *  o_  *  r  .  A(u,v)  then, 

A  D  L 


s  »  I  {(A  cos8 .  +  B  sin8.)  +  j(B  cos0.  -  A  sin0.)} 
m  i  i  J  i  i 

u,v 

The  computer  program  listing  which  incorporates  this  implementation 
of  random  reflectivity  is  presented  in  Appendix  7.1.  The  listing 
includes  subroutines  for  random  number  generation  which  were  supplied  by 
the  University.  The  function  NRHRAN  generates  a  normally  distributed, 
uncorrelated  random  number  of  specified  mean  and  variance  equal  1.0. 
Zero  mean  was  chosen,  and  the  standard  deviation  is  choosen  to  result  in 
the  desired  power.  The  reflectivity  is  correlated  from  pulse  to  pulse 
by  the  subroutine  SCATTER.  However,  the  reflectivity  of  each  resolution 
cell  is  made  independent. 


'  '  rre :  j:  ed  Random  Reflectivity 
•;  g’ire  * .  5 


.4.  i  s;mi  i.a:mn  hesm.ts 

Several  different,  targets  were  used  for  ISAR  image  generation  to 
observe  the  effects  of  correlated  random  reflectivity  and  additive 
noise.  First,  an  image  was  generated  from  a  point  target,  Located  in 
tiie  center  of  the  field  of  view,  using  a  correlation  interval  length  of 
-4.0  puise  intervals,  and  no  added  noise.  This  image  is  on  the  left  (at 
i:i  Figure  -4.3.  A  correlation  interval  of  A.O  means  a  correlation  time 
. / .  :  .  ’  :  -•  p  •  is  .  I  *  ;  -■  r  h <\ t  r  in-*  r • '  i  nr  f  -i r'^^t  :  s  rt* <.  ■  u i ~ 

ible  line  to  t he  random  reflectivity.  Also,  the  .mage  noise  is  concon- 

he?  ween  bursts.  So,  the  image  is  less  random  in  the  delay  a i mens  ion , 
-:i :  -h  :  •••  'he-  no  r  i  Oi;r.  t  a  ;  ir.  Figure  .  j .  Fhe  image  mi  "no  r :  g:tr 
generated  from  the  same  point,  target  using  the  sane  -o rre; at  Lnterva  i 

:  :;c  .  uo  .  :ig  i  id  i  t  :  ve  ;v  i  se  using  an  SNR  el  10  dec  i  L>e  Is  .  it  :  s  -een  •  a' 


r no  added  noise  has  little  effect  on  the  image. 


For  a  comparison,  an  ISAR  image  was  generated  from  the  spherical 
target,  using  correlated  random  reflectivity  of  correlation  interval  0.5 
pulse  intervals,  and  additive  noise  with  SNR  of  0  decibels.  This  image 
is  shown  on  the  left  (a)  in  Figure  A. 4.  On  the  right  (b)  in  the  same 
picture  is  an  image  resulting  from  processing  a  mixer  signal  containing 
only  additive  noise.  The  result  is  completely  random  noise  ana  is 
*'  *  <?  s  ^  i  r  su  i  ^  r  »i  b  *  ^  f  r ■ '  m  r*.t?  •? v h <?  v <  i  h ■”} c* o 

Several  images  were  generated  in  a  parametric  study  of  the  effects 
of  the  correlation  interval  value.  A  point  target  in  the  center  :f  ‘he 
image  as  used.  No  added  noise  was  put  in  the  mixer  output  signal.  Six 
ime  were  generated,  each  with  a  (different  value  of  correlation 


S' 


?.VI 


‘'vl 


interval.  The  results  are  displayed  in  Figure  A. 5  with  horizontal  being 
the  delay  dimension.  The  correlation  interval  value  increases  from 
image  to  image  going  from  top  left  (a),  to  top  right  (c),  to  lower  left 
(d),  to  lower  right  (f).  The  noise  in  the  image  seems  to  progressively 
narrow  horizontally  into  a  vertical  line,  and  then  the  line  shortens 
toward  a  point,  as  correlation  interval  increases.  The  six  correlation 
interval  (Cl)  values  are  (a)  Cl  =  0.5,  (b)  Cl  =  4,  (c)  Cl  =  10,  (d)  Cl  = 
128,  (e)  Cl  =  256,  and  (f)  Cl  =  1100.  The  noisy  image  becomes  a  line 
when  the  correlation  interval  equals  128  pulse  intervals,  the  same 
length  as  a  burst,  the  point  at  which  reflectivity  is  essentially  random 
only  from  burst  to  burst. 


ISAK  Images  of  Point  Target  Using 
Various  Correlation  Interval  Values 
Figure  4.5 


v  *.*  *  *  ■* 
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The  simulations  presented  in  Figure  4.5  were  repeated  using  the 
sphere  described  previously  as  the  target  to  be  imaged.  The  results  are 
presented  in  Figure  4.6.  It  is  seen  that  as  correlation  interval 
increases  the  target  becomes  more  recognizable. 


(a) 


(b) 


(C) 
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Additionally,  examples  of  windowed  ISAR  images  were  generated  to 
see  the  effects  of  windowing  on  diffuse  targets  and  additive  noise.  On 
the  left  (a)  in  Figure  4.7  is  a  windowed  ISAR  image  of  a  center  point 
target  having  a  correlation  interval  value  equal  1100  pulse  intervals 
and  no  additive  noise.  Windowing  seems  to  make  this  point  target  look 
like  an  amorphous  blob.  The  image  on  the  right  (b)  in  the  same  figure 
has  the  same  target  with  additive  noise  of  SNR  equal  0.0  decibels  added 
in  the  signal.  The  resulting  image  has  additional  blurry  background. 
These  images  suggest  that  using  windowing  in  generating  ISAR  images  in 
the  presence  of  noise  does  not  improve  image  quality. 


(a)  (b) 

Windowed  ISAR  Images  of  Center 
Point  Target 


Figure  4.7 
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4.4  POWER  SPECTRUM  ESTIMATION 

1 

As  described  before,  the  mixer  output  from  a  diffuse  point  target, 

with  no  added  noise,  is  a  temporally  correlated  random  reflectivity 

I 

multiplied  by  a  delay-doppler  signal.  In  generating  an  image  from  this 

P 

mixer  output,  Fourier  transforms  are  calculated.  The  result  can  be 

1 

1 

viewed  as  the  convolution  of  the  power  density  of  the  target,  the 

1 

desired  image,  with  the  Fourier  transform  of  the  random  reflectivity. 

1 

The  transform  squared  magnitude  of  a  sequence  of  random  numbers  is 

called  a  periodogram  in  (6).  Also,  this  is  described  as  an  estimate  of 

H 

the  power  spectrum  of  the  random  process.  If  the  Fourier  transform  were 

! 

P 

an  unbiased  and  consistent  estimate  of  the  power  spectrum,  the  error  in 

| 

the  estimate  of  the  power  spectrum  would  approach  zero  as  the  number  of 

1 

samples  of  the  random  process  gets  very  lar^e-  However,  this  is  not  the 

R 

case.  It  is  shown  in  (6)  that,  the  periodogram  is  a  biased  and  incon- 

P 

sistent  estimate  of  the  power  spectrum  of  a  random  process.  The  esti-  1 

K 

equal  the  power  spectrum  using  limited  sample  length,  and  the  estimate 

w 

is  inconsistent  because  the  variance  of  the  periodogram  does  not 

§ 

approach  zero  as  more  samples  are  used  in  the  calculation.  This  means 

I 

that  1SAR  processing  as  it  is  conventionally  performed  cannot  accurately 

reconstruct  the  power  density  of  of  a  diffuse  target,  equivalently  the 

i 

scattering  function.  An  approach  for  improvement  of  the  power  spectrum 

i 

estimation  is  the  Welch  method  described  in  (6).  In  this  method,  the 

V 

data  is  windowed  before  the  periodogram  is  computed  and  several  of  these 

1 

smoothed  periodograms  are  averaged.  This  would  require  the  averaging  of 

1 

several  reconstructed  images.  The  Welch  method  yields  a  power  spectrum 

msmmm 

V 


estimate  which  is  biased  for  limited  sample  length,  though  asymptoti¬ 
cally  unbiased,  and  the  variance  of  the  estimate  is  expressed  as 
follows . 

Variance  [B  (w)]  ~  1^  P^(w)  , 

K 

where  B  is  the  Welch  estimate,  K  is  the  number  of  periodograms  averaged, 
and  P  is  the  power  spectrum  of  the  random  process.  The  Welch  method 
sacrifices  spectral  resolution  and  bias  for  a  consistent  spectral 
estimate.  Using  the  Welch  method  with  limited  data  probably  will  not 
result  in  a  good  reconstruction  of  a  scattering  function. 


5.  SUMMARY  AND  CONCLUSIONS 


This  thesis  work  consisted  of  developing  an  ISAR  imaging  simulation 
using  a  computer  program  to  input  a  target  scattering  function,  generate 
a  simulated  mixer  output  signal,  and  reconstruct  an  image  of  the  target 
from  the  mixer  output  signal.  In  the  results  of  Chapter  3,  the 
reflectivity  process  was  fixed  at  a  constant  value  for  each  resolution 
cell,  allowing  the  verification  of  the  ISAR  processing  since  the 
reconstructed  image  is  expected  to  look  like  the  scattering  function  of 
the  target. 

With  the  ISAR  simulation  verified,  the  next  step  reported  in 
Chapter  4  was  to  implement  random  reflectivity  and  additive  noise.  This 
was  to  use  the  slowly  fluctuating  target  model  described  in  (5).  The 
complex  random  reflectivity  process  was  correlated  in  time  using  a 
recursive  filter.  The  random  reflectivity  was  kept  spatially 
independent.  The  results  of  generating  ISAR  images  using  random  reflec¬ 
tivity  of  diffuse  targets  and  additive  noise  are  presented.  The 
reconstructed  images  show  that  the  random  reflectivity  has  great  effect 
on  the  image.  As  a  result,  the  actual  target  cannot  be  recognized  from 
the  image  unless  the  random  reflectivity  is  highly  correlated. 

The  conclusion  is  that  ISAR  processing  attempts  to  measure  the 
power  density  of  the  mixer  output  signal  which  would  look  the  same  as 
the  scattering  function,  but  the  random  reflectivity  makes  the 
calculated  power  density  a  random  quantity  and  so  the  result  is  a  noisy 
image.  Fourier  transform  techniques  yield  unsatisfactory  results  in 
this  application  of  power  spectrum  estimation  theory. 
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APPENDIX  7.1 
Program  Listings 

The  listings  of  the  Fortran-77  programs  and  subroutines  used  to 
simulate  the  mixer  output  signal  and  perform  ISAR  processing  to 
reconstruct  an  image  are  included  in  this  appendix.  A  brief  description 
of  each  program  and  subroutine  follows. 

PR0C3,  listed  in  Figure  7.1,  calculates  a  mixer  output  signal  from 
an  input  scattering  function  data  file.  A  non-statistical  target  is 
assumed.  Subroutine  TARGET  reads  in  the  scattering  function. 

VT  *  Tangential  velocity  of  the  target 
XFOV,  YFOV  »  Image  field  of  view  dimensions 

CRX,  CRY  *  location  of  center  of  rotation  in  scattering  function 
T1  *  pulse  width  (sec) 

T3  *  Burst  length  (sec) 

W  *  angular  rotation  rate  in  radions  per  second 
BW  ■  bandwidth  of  burst 
FO  *  lowest  frequency  (Hz) 

FC  ■  center  frequency 
DF  *  frequency  step  size 
WL  »  wavelength 

N  “  number  of  pulses  and  bursts 

XRES  *  resolution  in  meters 

RWCF  *  range  walk  correction  factor 

Subroutine  IMAGE2D,  listed  in  Figure  7.2,  performs  ISAR  processing 
to  reconstruct  an  image  from  radar  received  signals.  The  discrete 
Fourier  transforms  are  implemented  using  Singleton's  Fast  Fourier 
transform  algorithm. 
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PR0C4,  listed  in  Figure  7.3,  calculates  a  radar  received  signal 
from  an  input  scattering  function  data  file.  This  routine  is  a 
modification  of  PR0C3  to  include  the  effects  of  a  statistical  target  and 
additive  noise. 

Cl  *  correlation  interval  (#  samples)  of  random  reflectivity 
SNRCfi  *  signal  to  noise  ratio  for  additive  noise  in  decibels 
Subroutine  SCATTER,  listed  in  Figure  7.4,  generates  components  of 
complex  random  reflectivity  which  are  correlated  in  time  by  recursive 
filtering  white  gaussian  noise. 

Subroutine  IMAGE2DW,  listed  in  Figure  7.5,  uses  data  "windowing"  as 
part  of  the  Fourier  transform  techniques  which  generate  an  image. 


55- 


PRCGRAM  PROC3 
DIMENSION  S'  1 28 1  , O'  1 23) 

REAL  LAMBDA! 128. 128 ) 
Coop  I f  x  FRAME  <236,25*) 
Int*g»r«4  tun* 

:  wmon  btKI.  'uns 


C» 

C*  ISAR  IMAGING  USING  A  STEPPED-FREQUENCY  WAVEFORM 

C» 

C»  BOB  LEWIS,  DECEMBER  1986 

C *••••••••••••••••••••••••••••••••••••• ••••••••••••••*• 

03.  E8 
C2-2./C 
PI-3.1413927 
pi  2*2  .  *P! 

RES-0 .o 

N-12B 

Y RES- RES 

XRES— RES 

WLCF-0 . t3 

FC-C/ULCF 

RWC— PI 2/YRES/N 

XFOU-N-XRES 

YFOU— N*YRES 

CRX— <N/2-l ) *XRES 

CRY-<N/2-l >  *YRES 

Tl  — 2  .  *YFOV-'C 

TI2-TI/2 

8WC/<2.«YRES> 

F0-FC-BU/2 

WL-C/F9 

TMETA-UL/! 2*XRES) 

RANGE- 1  2. 5* 1 832. 

UT-400.*1832./3d0e. 

W-VT/RANGE 
T3— THETA/ <  W*N ) 

T2-T3/ <N-1 . 

T4-N*T3 
OF— BW/ < N- 1 ) 

CALL  i cr **t < ' En t*r  output  fil*n*m*  '.lun*> 


REAO  IN  THE  TARGET' S  SCATTERING  FUNCTION  •**•** 
CALL  t»pg*t(  LAMBDA) 

c******  CALCULATE  MIXER  OUTPUT  SIGfsALS  AND  STORE  IN  FRAME 
L— »5 

DO  180  T-e . ,T4,T3 
DO  30  K-t ,128 
S' Ki-9 , 

Q<  K)-0 . 


30  CONTINUE 

DO  30  J-l .128 

<m<  .J-l  i  •XRES-CR* 

*  .  »u 


PR0C3 

Figure  7.1 


VC 2 -U«C2 
«UCF-U»T*RUC 
00  70  I«l ,128 

if<lambda< i , j> .eq.0 . >go  to  ?e 

A- LAMBDA < 1 , J)*«0 .3 
Y-( 1-1 >*YRES 
DELAY-C2-Y 
DO  60  K-l ,128 
Kl-K-1 
FI-K1*DF*F0 

PH I -PI 2«FI •<  DELAY -VC2*T»  »K1 *RUCF 
S<K)-S(K)*A*COS<PHI ) 

Q( K)-0( K ) »A«SIN< PHI > 

CONTINUE 
CONTINUE 
CONTINUE 
DO  90  K— l,12B 

FR**IE<K*64.L>-CMPLy<SCK)  ,Q(K)  ) 

CONTINUE 
L-L*l 
100  CONTINUE 

C»>«**«  NORMAL  1 2E  MIXER  OUTPUT  SIGNAL  TO  PEAK- 1  . 0 
PX-0  . 

DO  120  L-1,236 
DO  110  K-l .236 

PK-MAX  <  PK  ,  ABS  (  FRAME  C  L  ,  K  M  > 

110  CONTINUE 

120  CONTINUE 

PKINV-l  ,/PK 
DO  140  L-1,236 
DO  130  K-l ,236 

FRAME r  L  ,  K  '  —FRAME  (  L  ,  K  1  •  PKINN1 
130  CONTINUE 

140  CONTINUE 

C»**»*«  RECONSTRUCT  THE  TARGET'S  IMAGE  •****•**•• 
CALL  i m*g*2D<  FRAME ' 

•  top 
*nd 


60 

70 

80 


90 


PR0C3 

Figure  7.1  (Continued) 
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Sufcroutio*  .  mag»2C"  •  r  am*  ; 

R*al  . mag* (236, 236) 

Camp  1 *x  f ( 234) 

Complex  f run* <236 ,236) 

common  'plkl  '  1yn» 


This  Subroutine  reconstruct*  t  radar  image 
by  uting  Fourier  transform  technique* 

Bob  l*m>s.  December  if 86 
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r*8 

•«  DELAY  DIMENSION  RANGE  PROFILES  *•« 
00  120  j-1.236 
DO  100  .*1.25« 

f < i  >*eonjg<  frame! i  ,j ) > 

CONTINUE 
CALL  FFT(f.M) 

DO  1  10  i*l  .236 

f ram# <  i  .  j  >*conjg( f < i > ) *0 .8037063 
CONTINUE 
CONTINUE 


>•  DOPPLER  DIMENSION  PROFILES 
00  138  1*1.236 
00  123  j-1,236 

f<j)*cmplx<0. ,0.) 

CONTINUE 
DO  130  J-l .64 

f < j»o4)*fram*( i , j«128) 
f < ,rM28>*fr*m*< i , j +64 ) 
CONTINUE 
CALL  fft(f.M) 

00  140  j-t ,128 

imag*< i  . j*128)»»b*( f (j ) >**2 
. mag* <  . , j >*ab»'. f ( j ♦ 1 28' < **2 
CONTINUE 
CONTINUE 


nby t*s* 236*236*4 

CALL  fwr  i  t*(  I  un*  ,  imag*  ,  nbxtes) 

R*  turn 

END 
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PROORAH  PR0C4 
REAL  LAM«OA< 128,128) 

Int#Q*r*4  1  RANOl  ,  I  RAN02 

Como i**  FRAMED  236. 236 > . CNO IS 

Common  o  1  K  I  >•  I  uni 

Common  /OIK 3/  Y 1 , Y2 , I RAND  I  , 1 RAND2 
Oim«nsion  GAUSS ( 2) 

C •••*•••***•••***••*•*•••*********•••«•«•*•****• •*****•• 
c* 

C»  I  SAP  WITH  CORRELATED  COMPLEX  RANDOM  REFLECTIVITY 

C*  ANO  ADDITIVE  COMPLEX  GAUSSIAN  WHITE  NOISE 

C*  SOS  LBWIS.0CC8M8ER  I  *84 

C-3.E9 

C2-2./C 

PI-3.1413P27 

PI  2- <4. 283 1834 

RES-0.6 

YRES-RES 

XRES-RES 

SMI  28 

ULM  .83 

FC-C/UL 

XFOV-m«xRES 

YFOV-N«YRES 

CRX*<N/2-l >«XRES 

CRY*<N/2-l >»YRES 

Tl*2.«YFOV/C 

T12-T1/2. 

*MC/<2.»YRES> 

Fi-FC-SU/2 

THrrA*2MS!N(WL/(4*XRES>  ) 

RANGE* 1 2. 3*1 892 . 

VT*488.«1892./348t . 

LMVT/RANGE 

T3*THETA/<W*N) 

T2-T3/<N-1 ) 

T4*N«T3 
OF-BU/(N-l ) 

PRINT  •, 'Input  corrtlition  intrvtl 
REAO  «,CI 

CALL  f cr*»t < ' En t*r  output  ,lunt) 

C*»***«  REAO  IN  THE  TARGET’S  SCATTERING  FUNCTION 
CALL  t  «r  g«  t  <  LAM80A ) 


SNRO^«.8 
I  RAND I *3224?  l  21 * 
I  RAND 2-7 


PROC4 


Figure  7.3 
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sun  TOTAL  POWER  IN  SCATTERING  PNC  LAMS  DA 
INITIALIZE  20  ASSAY  CONTAINING  MIXER  OUTPUT  SIGNALS 


00  4*  1-1  ,128 
DO  38  J-l  ,  1  28 

SUMP" SUMP *LAMGOA! I , J) 
PSAMC < I ♦ 44 , J* 44 >•* . 

3#  CONTINUE 
M  CONTINUE 


•••*••  CALCULATE  MIXER  OUTPUT  SI94ALS  ANO  STORE  IN  FRAME  ••••••• 

Loops  PS  k  IBS  change  resolution  cell 
■••*•••  Loop  9f  steps  burst  R  ,  Loop  AS  stsps  puls*  R  ••••••• 

’*••••*  Subroutine  SCATTER  picks  rindor  reflectivity  •••**•• 


DO  1SS  >1,128 

X"! J-l > sXRES-CRX 
U"XeWeC2 
DO  PS  1-1,128 

I  F I LAMBDA! I , J ) . EO . I . ) GO  TO  PS 

G"!  LAMBDA ( I , J)e* .3)*M .5 

Y*< I -1 )*YRES 

DELAY-C2*Y 

L-43 

CALL  NSMRAN ( I  RAND 1 , I  RAND 2 , GAUSS  > 

Y1 "GAUSS! 1 >*G 
Y 2-GAUSS ! 2>  *G 
DO  8S  T-p. ,T4,T3 
DO  AS  K-l ,128 

F1-<K-1 >#DF*F8 
PHI*P12«FI*CDELAY-W«T> 

CALL  SCATTER! 3, Cl ,A,B> 

S-A#COS!PMI )-8*SIN!PHI > 

>B*COS!  PHI  >*A«SIN<PMI  > 

FRAME ! L , K»  A4 ) “FRAME  I  L .  X  ♦ AA  > ♦ CMPLX ! S , 0 ) 
CONTINUE 
L»L»  l 
CONTINUE 
CONTINUE 


•****•  Calculate  standard  deviation  o4  additive  noise  Aram  SNR 
••••••  Pick  additive  white  Gaussian  noise  values 

•*••••  notas! I te  si xer  output  si gnal iFRAME)  to  power-i .1 

•*•••*  Add  noise  to  mixer  output  signal 


ACS- 1 ./SUMPeef .3 

lS.eelSNROS^I*. > 

ANO  IS- 1 ,/SNReeS . 9 
SDNO I  S— ANO  I  S«S  .  7871 
DO  l 28  L-A5.IP2 
DO  US  K-A3.IP2 

CALL  NRMRAN!  IRANOl  ,  IRAN02, GAUSS) 
SNO 1 SmGAUSS ! I ) *SONO I S 
a«l  S-GAUSS I  2  >  eSONO  I S 
CNOI S-CMPLX !  SNO  I S ,  GNOI S  > 

FRAME < L , K > -PRAME ! L , K ) RAC S » CNO I S 
US  CONTINUE 
1 28  C94TINUE 


CALL  image 20! FRAME) 
STOP 

ENO 
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Figure  7.3  (Continued) 
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Ceeeeee  Thi*  uibrsutmt  generate*  complex  random  reflectivity 
£*••**•  which  are  correlated  in  time  by  reeureive 

;***•*•  *  > ' ter i ng  white  Gao  »» .  an  no i *e ,  »upp I i ed  B>  NPMfiAN . 
*••••••  Cl *c or r e l at i on  interval  ie  tamp  let; 

Ceeeeee  Bob  Lewie,  Dec.  If** 


Integer**  i randl , i rand2 

common  /blK3/  yt ,y2, i randl , i rand. 

dimeneion  gau**(2> 

CALL  nrmran< i randl  , i rand2 ,gauee> 

* l“gaues< 1 5  *g 

al pha"-l .  t/e i 
rho»exp  < alpha) 
gamma“( 1 ,-rho**2>**8 .3 

y 1 “r ho*y 1 « gamma *x 1 
y2“rho*y2»gamma*x2 


re  turn 
ENO 


SCATTER 
Figure  7.4 
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Subroutme  image20<  f  ram*  ) 
»»*'  i mag* ( 236 , 23* > 

Cofttt  i  *  *  *  1  2*6  > 

Complex  * r erne (23*, 236 > 
dimension  wl(128) 
common  /plKI/  I unt 


C< 

C<  This  Subroutine  reconstructs  *  r»d*r  imag* 


C»  by  using  mi n flowed  Four i er  transforms 

C<  Sob  Lewis,  December  178* 

CX<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 

H-8 

sf-3.761 

CONSTRUCT  WINOOU  #**•*•*•*•*•***••*( 

00  10  i «1 .128 

wl ■ . »f<*«pi-0.09<< i -o4 • e*2 • 

18  CONTINUE 

Ceeeeeee  DELAY  DIMENSION  RANGE  PROFILES  •  < 
DO  128  j-1,236 
00  188  i-l,236 

f < i >-con j  g<  4 r am*  < i  , j  > ’ 

108  CONTINUE 

DO  181  1-1,128 

4 < i ♦ 64)— f  < i *64)  ewl (  <  > 

181  CONTINUE 

CALL  FFT <  4 ,M> 

DO  118  i-l ,236 

4 r ame< I  , j  >-Conj g<  4 < I  > ) <8.8837863 
110  CONTINUE 

128  CONTINUE 

C«««<  OOPPLER  DIMENSION  PROFILES  «* 

00  138  i«t,236 
00  123  u-1 ,236 

4 <  j >-cmp 1 x  <  8  .  ,  8  .  > 

123  CONTINUE 

DO  138  u-1 ,64 

4 <  j<64>— frame < i  , j* 1 28) 
f < j< l 28>-f r am* c i ,j<64> 

138  CONTINUE 

00  131  j-l ,128 

4 <  j»64>-f < J»64) ewl <  j  > 

131  CONTINUE 

CALL  fft(f,M> 

DO  140  j-l ,128 

image ( i , j ♦ 1 28>-«bs( f  <  j ) > <<2 
i mage < i ,j>— *b*(f<j<128>>«2 
140  CONTINUE 

130  CONTINUE 


nDy tes— 236*236<4 

CALL  fwr i te< 1 uns . i mage , nby tes > 
Return 

END 
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Figure  7.5 
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ABSTRACT 

The  expectation-maximization  algorithm  for  computing  maximum-likelihood  estimates 
iteratively  is  used  to  develop  a  new  approach  for  processing  inverse  synthetic-aperture  radar 
data  to  form  images  of  fluctuating,  diffuse  radar-targets.  The  scattering  function  of  the  target 
is  imaged  by  jointly  estimating  the  power  spectra  of  wide-sense  stationary  reflectivity-processes 
occurring  in  all  the  range  cells  that  span  the  target.  The  complex-valued  reflectivity  processes 
are  also  estimated.  The  results  we  obtain  apply  to  imaging- radar  systems  operating  at  radio  and 
optical  frequencies  when  target  echos  have  no  specular  or  glint  components. 


t  This  work  was  supported  by  the  Office  of  Naval  Research  under  contract  N00014-86-K0370. 


1.  Introduction 


An  inverse  synthetic-aperture  radar  (ISAR)  system  is  used  to  form  an  image  of  a  radar 
target.  This  is  accomplished  by  illuminating  the  target  with  a  series  of  pulses  and  observing  the 
return  echos.  Each  patch  on  the  target  introduces  a  certain  amount  of  propagation  delay  and 
doppler  shift  to  a  pulse  it  reflects,  the  amount  depending  on  the  range  and  range  rate  of  the 
patch  relative  to  the  radar  system.  The  received  signal  for  each  illumination  is  a  complicated 
superposition  of  the  reflections  from  all  patches  that  make  up  the  extended  surface  of  the 
target.  The  goal  in  processing  the  received  signal  is  to  produce  an  image  of  the  target. 

The  design  of  an  ISAR  system  includes  the  selection  of  the  transmitted  waveform,  the 
selection  of  a  model  for  the  reflection  process  by  which  a  portion  of  the  transmitted  waveform 
is  returned  to  the  receiver,  and  the  selection  of  the  algorithm  used  to  process  the  received 
waveform  to  create  the  target’s  image.  The  beam-width  of  the  radar  antenna  relative  to  the  size 
of  the  target  is  also  a  design  consideration;  images  can  be  produced  by  either  scanning  a 
narrowly  focused  beam  over  the  target  in  some  type  of  raster  pattern  or  by  illuminating  the 
entire  target  in  spotlight  mode  with  a  wide,  relatively  unfocused  beam.  Our  concern  will  be 
with  forming  images  of  rotating,  rough  targets  using  a  spotlight- mode  radar. 

Stepped- frequency  and  l  inear- FM  chirp  are  two  modulation  formats  used  with  transmitted 
pulse-sequences  in  spotlight-mode  radar- imaging.  Stepped-frequency  pulses  are  described  by 
Prickett,  Wehner,  and  Chen  [1],  Ruttenberg  and  Chanzit  [2],  and  Chen  and  Andrews  [3].  The 
target  is  illuminated  with  a  sequence  of  N  pulse-groups,  where  each  group  is  identical  and  is 
formed  from  a  sequence  of  N  narrow  pulses  having  equal  incremental  steps  in  frequency.  The 
usual  approach  for  processing  delay-  and  doppler-shifted  echos  acquired  by  illuminating  a 
rotating  target  with  this  waveform  consists  of  two  steps.  The  first  is  to  place  the  data, 
consisting  of  one  sample-value  per  transmitted  pulse,  into  delay  (or,  range)  bins  by  separately 
Fourier  transforming  the  N  sample  values  from  each  pulse  group.  The  resulting  delay -binned 
data  are  placed  in  the  rows  of  an  NxN  matrix,  where  each  row  contains  the  transformed  data 
from  one  pulse-group.  In  the  second  step,  the  columns  of  this  matrix  are  Fourier  transformed 
to  obtain  a  doppler  (or,  cross- range)  profile  at  each  delay.  The  resulting  two-dimensional  array 
is  intended  to  be  the  target’s  complex-valued  reflectance  function  in  delay  (range)  and  doppler 
(cross-range)  coordinates,  the  magnitude  or  squared  magnitude  of  which  can  be  displayed  as  the 


target’s  image.  This  processing  based  on  two-dimensional  Fourier  transforms  is  derived  using  a 
deterministic  analysis  that  does  not  account  for  statistical  properties  of  the  reflectivity  or  any 
noise  that  may  be  present. 

Wideband  "chirp"  pulses,  having  an  instantaneous  frequency  that  varies  linearly  with  time, 
are  also  used  for  radar  imaging.  The  common  approach  is  to  transmit  a  series  of  pulses,  each  of 
which  has  an  identical  envelope  and  chirp  rate.  A  variety  of  processing  approaches  are 
described  in  the  literature,  including  two-dimensional  Fourier  transformation  by  Mensa  [4]  and 
Walker  [5]  and  tomographic  reconstruction  by  Mensa  [4]  and  Munson,  O’Brian,  and  Jenkins  [6]. 

A  variant  of  the  stepped-frequency  format  in  which  each  narrow  pulse  in  a  pulse  group  is  a 
chirp  is  described  by  Blahut  [7] 

Bernfeld  [8]  introduced  the  concept  of  using  chirp- rate  modulation  with  processing  based 
on  tomographic  reconstructions  to  image  the  target's  scattering  function  for  radar  signals  having 
a  large  time-bandwidth  product.  With  his  approach,  the  target  is  illuminated  by  a  sequence  of 
linear-FM  chirped  pulses,  with  each  pulse  having  a  distinct  chirp  rate.  Bernfeld  notes  that  the 
output  of  a  matched- filter  receiver  for  a  radar  waveform  with  an  infinite  time-bandwidth 
product  is  a  line  integral  through  the  scattering  function,  where  the  angular  orientation  of  the 
line  of  integration  for  a  pulse  is  a  function  of  the  chirp  rate  of  that  pulse.  This  observation 
suggests  using  the  same  algorithms  as  used  in  x-ray  tomography  to  determine  the  scattering 
function.  Snyder,  Whitehouse,  Wohlschlaeger,  and  Lewis  [9]  extended  Bernfeld’s  approach  to 
include  waveforms  with  more  modest  time- bandwidth  products.  This  is  accomplished  by  noting 
similarities  to  the  tomographic  imaging  of  radioactive  tracers,  where  ideal  line-integrals  are  not 
available. 

The  approach  we  describe  in  this  paper  can  accommodate  the  stepped-frequency  and  chirp 
formats  as  well  as  others.  It  represents  a  continuation  of  our  examination  of  how  the 
approaches  currently  used  in  emission  tomography  can  be  applied  to  radar  imaging.  The 
method  we  have  described  in  [9]  is  adapted  from  the  best  algorithm  for  time-of-flight  emission 
tomography  when  the  processing  of  the  data  is  required  to  be  linear  [10].  A  more  fundamental 
approach  described  by  Snyder  and  Politte  [11]  leads  to  nonlinear  processing  and  improved 
accuracy  in  forming  images  of  radioactivity  distributions.  In  this  paper,  we  describe  the  initial 


results  we  have  obtained  in  adopting  an  analogous  approach  for  radar  imaging.  This  relies  on 
the  use  of  a  mathematical  model  for  data  acquired  in  a  spotlight- mode  radar  and  on  the  use  of 
the  method  of  maximum-likelihood  estimation. 

Van  Trees  [12]  and  Shapiro,  Capron,  and  Harney  [13]  describe  models  for  fluctuating, 
diffuse  radar-targets.  The  models  are  based  on  the  assumption  that  the  surface  of  the  target  is 
rough  compared  to  the  wavelength  of  the  radar  illumination.  The  reflectance  is  modeled  as  a 
Gaussian  wide-sense  stationary,  uncorrelated  scatter  (WSSUS)  process,  which  is  uncorrelated  in 
range  and  temporally  stationary.  Such  models  are  accurate  for  some  targets  in  radio-  and 
optical- frequency  radar-imaging  systems,  but  not  all.  An  important  effect  not  included  in  our 
present  results  is  that  of  glint  or  specular  components  in  radar  returns.  We  arc  currently 
attempting  to  extend  our  approach  to  include  these  additional  effects. 

Model  based  approaches  that  use  statistical  estimation-theory  techniques  appear  less 
frequently  in  the  large  liteiature  of  radar  imaging.  One  example  is  that  of  Frost,  Stiles, 
Shanmugan,  and  Holtzman  [14],  who  use  a  multiplicative  model  and  Wiener-filtering  techniques. 
Our  approach  differs  in  that  the  model  we  adopt  for  the  return  signal  is  more  complicated  than 
a  simple  multiplicative  one  and  depends  explicitly  on  the  transmitted  waveform  through  a 
spatial  integration  over  the  target.  We  also  do  not  restrict  the  processing  to  be  linear,  in 
particular,  we  show  that  the  processing  of  the  received  data  for  producing  the  maximum-likeli¬ 
hood  estimate  of  the  target’s  scattering  function  is  not  linear. 

Radar  imaging  systems  generally  produce  estimates  of  one  or  the  other  of  two  quantities 
that  can  be  viewed  as  the  target’s  image.  Some  approaches  produce  an  estimate  of  the  target’s 
complex- valued  reflectivity  as  a  function  of  range  and  cross  range,  while  others  produce  an 
estimate  of  the  target’s  scattering  function.  For  the  new  approach  we  describe,  the  reflectivity 
process  is  modeled  as  a  complex- valued  Gaussian  random- process  that  is  temporally  stationary 
and  spatially  white.  The  scattering  function  is  the  power-density  spectrum  of  this  process  as  a 
function  of  delay.  We  treat  both  the  reflectivity  process  and  its  second-order  statistic,  the 
scattering  function,  as  unknown  quantities.  The  iterative  approach  we  develop  yields  the 
maximum-likelihood  estimate  of  the  scattering  function  and,  simultaneously,  the  condition¬ 
al-mean  estimate  of  the  reflectivity  based  on  statistics  which  are  consistent  with  the  estimated 


scattering-function.  Thus,  both  of  the  quantities  treated  separately  in  other  imaging  schemes 
are  produced  simultaneously  with  our  new  approach.  This  is  unique  to  our  approach,  and  we 
feel  that  it  is  important. 

We  will  develop  a  necessary  condition,  called  the  trace  condition,  which  the  maximum-like¬ 
lihood  estimate  of  the  target’s  scattering  function  must  satisfy.  This  equation  appears  to  be  very 
hard  to  solve  directly.  As  a  consequence,  we  reformulate  the  imaging  problem  using  the 
concept  of  incomplete-complete  data  spaces  and  then  use  the  expectation-maximization 
algorithm  to  derive  an  iterative  algorithm  for  producing  the  maximum-likelihood  estimate  of  the 
scattering  function.  This  procedure  also  yields  the  conditional-mean  estimate  of  the  reflectance. 
The  technique  we  use  to  accomplish  this  parallels  that  described  by  Miller  and  Snyder  in  [15] 
for  power-spectrum  estimation  and  extends  their  results  to  include  indirect  measurements  of  the 
process  whose  spectrum  is  desired;  the  process  is  now  measured  following  a  linear 
transformation  and  in  additive  noise.  As  shown  by  Turmon  and  Miller  [16],  this  approach  to 
spectrum  estimation  results  in  estimated  spectra  with  less  bias  and  mean-square  error  than  other 
approaches  discussed  in  the  literature.  We  expect  similar  improvements  will  be  seen  in  radar 
imaging  of  scintillating,  diffuse  targets  when  this  new  technique  is  used. 


2.  Model 


The  model  and  notation  we  shall  use  closely  follows  that  of  Van  Trees  [12,  Ch.  13].  The 
complex  envelope  of  the  transmitted  signal  will  be  denoted  by  2£t1/jst(0  where  £T  is  the 
transmitted  energy  and  s^(t)  is  normalized  to  unit  energy.  Thus,  the  transmitted  signal  is  given 
by 


yfTE  TRe[sT[t)exp{j2nf  0t)], 


(1) 


1 

1 


where  /o  is  the  carrier  frequency  of  the  radar.  For  a  stepped-frequency  waveform  consisting  of 
N  pulse  groups  with  N  pulses  in  each  group,  the  complex  envelope  is  given  by 


N-  I  A/-  1 

srU)«  Z  I  p(f-nrp-ir>xp(y2*d„(t-nrp-irj). 

1-0  -.-0  (2) 

where  Tp  and  Tt  denote  the  time  interval  between  pulses  in  a  group  and  between  groups, 

|  respectively,  A„  is  the  increment  to  the  carrier  frequency  fo  of  the  n -th  pulse  in  a  group,  and 

i 

p(t)  is  the  complex  envelope  of  a  pulse.  For  pulses  typically  assumed,  the  envelope  |/?(f)|  is  a 
narrow,  rectangular  function,  and  the  phase  arg[p(/)]  is  zero.  For  a  sequence  of  N  pulses  having 
chirp- rate  modulation, 

N- I 

srU)~  Z  p(<-iTp)exp(//r0,(r-iT,)2), 

‘■°  ’  (3) 


where  f)\  is  the  chirp  rate  of  the  i -th  pulse,  and  iTp  is  its  delay. 

Patches  on  the  target  with  a  two-way  delay  in  the  interval  [r,r+Ar)  reflect  a  signal  that  is 
incident  on  the  patch  at  time  /  with  strength  fe(/,r)Ar.  Consequently,  the  complex  envelope  of 
the  received  signal  srU)  following  the  illumination  of  the  target  by  sjO)  *s  given  by  the 
following  superposition  of  returns  from  reflecting  patches  at  all  the  two-way  delays  r 


U)  =  -J2 E T  J sT{t  —  r  )b(^t -  ,  r  jdr . 


(4) 


Van  Trees  [12]  and  Shapiro,  Capron,  and  Harney  [13]  discuss  the  reflectance  process  b{t.r) 
for  targets  that  are  rough  compared  to  the  wavelength  at  the  carrier  of  radio-  and 


optical-frequency  radar-systems.  For  such  diffuse  targets,  without  any  glint  components  in  the 
return  signal,  bit.r)  may  be  taken  to  be  a  complex-valued  Gaussian  random- process  with  zero 
mean  and  covariance 


E[b(t ,  r)6*(t ' ,  r  ')]  -  /C(t-  t',r)6(r-  r'),  (5) 

where  denotes  complex  conjugation,  the  impulse  function  in  delay  results  because  of  the 
assumption  of  uncorrelated  scattering,  and  K(tj)  is  the  covariance  function  of  the  reflectance  at 
each  delay  r.  The  Fourier  transform  of  K(t.r)  with  respect  to  t  is  the  scattering  function  S(/,r) 
of  the  target. 


K{t  ,r)exp(-  j2n  ft)dt . 


(6) 


This  is  the  power-density  spectrum  of  the  reflectance  process  for  all  scattering  patches  at  delay 

T. 


We  will  model  the  complex  envelope  of  the  total  return  signal  r{t)  as 

r(t)  -  s„(t)  u/(t).  (7) 

where  w(t)  is  complex- valued  Gaussian  white  noise  that  is  uncorrelated  with  the  reflectance 
process.  The  mean  of  w{t)  is  zero,  and  the  covariance  is 

E[u/U)t0*(t')]-Afofi(f-r).  (8) 

The  scattering  function  S(f,T)  of  a  diffuse,  rotating  radar  target  provides  an  image  of  the 
target  in  doppler  (or  cross  range)  and  delay  (or  range)  coordinates.  S(/,r)A/Ar  is  the 
mean-square  strength  (or  power)  of  the  reflectance  of  all  patches  on  the  target  having  a  doppler 
shift  in  the  interval  (/,/+A/)  and  a  delay  in  the  interval  [r,r+Ar).  We  may,  therefore,  state  the 
problem  of  imaging  a  diffuse  radar- target  as  that  of  estimating  the  scattering  function  S(/,t)  or, 
equivalently,  the  covariance  function  K(tj)  given  radar-return  data  {/•(/),  T,  <  t  <  T{ }  on  an 
observation  interval  (TiyT{). 


-  6 


discrete  model 


l 


to 


RS 
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In  anticipation  of  using  discrete-time  processing  of  radar  data  to  produce  images,  we  now 
state  the  discrete  version  of  our  model  as  follows.  We  are  given  N  samples  of  the 
complex- valued  radar-data  corresponding  to  (7), 


r(n)  -  s*(n)+  w(nj,  n-O.l . /V  -  1  , 

where  w(n)  is  a  white  Gaussian-sequence  with  zero  mean  and  covariance 

E[u/(/i)u/*(n')]-  N06„ 

and  where  the  signal  samples  corresponding  to  (4)  are  given  by 

s*(n)  -  Et  sT[n,i)b(n,  t ) ,  n-0,1 . N-  1 . 

In  this  expression,  we  define  si(n.i)  and  b(n,i)  in  terms  of  the  transmitted  signal  and  the 
reflectance  process,  respectively,  according  to: 


and 


b(rc,i)-b|  nAt  -  -i/3r  ,iAr  ]Jt, 


(9) 


GO) 


(11) 


(12) 


(13) 


where  A/  ard  Ar  are  the  sampling  intervals  adopted  in  the  discretization  in  time  and  delay, 
respectively.  We  assume  that  the  target  has  a  finite  extent,  so  b(n,i)  and,  therefore  also,  terms 

forming  the  sum  in  (11)  are  zero  for  i  outside  the  /  values  m,  m+1 . m+1-1  starting  from  the 

minimum  two-way  delay  corresponding  to  m.  This  discrete  reflectance  is  a  Gaussian  sequence 
with  zero  mean  and  covariance  given  by 


E[b(rt,ijb*(rt',i')]-K(n-rt',t)6((-. 

The  discrete  scattering  function  S(/,i)  is  the  Fourier  transform  of  K(n,i), 


(14) 


-  7  - 


(15) 


S(f,i)~  Y.  K{n,i)exp{-j2n/n). 

/l  ■  -  • 

The  imaging  problem  for  the  discrete  model  is  to  estimate  or  equivalently  the  covariance 

function  K(n.i),  for  all  frequencies  /  spanning  the  target  in  doppler,  and  for  all  delays  / 
spanning  the  target  in  delay,  given  the  radar  data  { r(n ),  n-0,  1 . N-I). 


matrix  model 

These  discrete  equations  may  conveniently  be  written  in  matrix  form  as  follows.  Define  r 
to  be  the  received-signal  vector  of  dimension  N, 


f  r(0)  \ 

1  r(l)  ' 


\r(N  -  1  )/ 


SK  +  W  ' 


(16) 


where  the  A^-dimensional  vectors,  sr  and  w,  are  given  by 


f  s*(  0)  \ 


s«(l. 


/  u/(0)  \ 

Ml) 


and  uj  ’ 


\sR(N-l)J 


(17) 


\w[N  -  1  )J 


Also,  define  5*  as  the  NIxN  rectangular  matrix  expressed  in  column-block  form  in  terms  of  / 
separate  NxN  matrices  according  to 


V-.J 


(18) 


where  Sj  is  an  NxN  diagonal  matrix  containing  sample  values  of  the  complex  envelope  of  the 
transmitted  signal  stO)> 


-  8 


/  s  r  ( 0 .  m  +  y ) 

0 

0  • 

0  \ 

0 

sr(  1  <rn*  j) 

0  • 

0 

0 

0 

• 

0 

V  0 

0 

st{N  -  1  ,m  +  j)J 

Further,  define  the  reflectance  vector  b  of  dimension  Nl  in  the  column-block  form  of  /  vectors 
according  to 


k 

k 


3? 


\b(/-  1 

where  each  b(i)  is  a  vector  of  dimension  N, 


b(0 ,  m  *  i) 
b(  1  ,  m  *  i) 


b[N  -  1  ,m*i) 


Using  (11)  and  these  definitions,  we  can  now  express  the  W-dimensional  signal  vector  5r  of  (15) 
and  (16)  as 

s h  -  j2ETS'b,  (22) 

where  "+"  denotes  an  Hermitian  transpose.  In  terms  of  these  defined  matrices,  the  received 
vector  has  zero  mean  and  covariance 


Kr  =  E(rr')  -E(s„s;)+E(u/u/*) 

=  2E TS~  t{bb')S  +  N  0 1 . 


I 


Then,  since 


E(b(t)b*(y))  =  K[i)6t  j, 


rrvm  jr« 


where  K(i)  is  the  Hermitian-symmetric  Toeplitz- matrix 


K  ( 0 ,  m  i )  K  *  ( 1 ,  m  + 1 ) 
K  ( 1  ,  m  +  i )  K[0  ,m  +  i) 

K  ( N  -  1  ,  m  *  i ) 


K*{N-  1  ,m  +  i 


K(0,  m  +  f ) 


it  follows  from  (23)  that  the  covariance  Kr  of  r  is  given  by 


Kr  -  2E  TS’ KS  +  A/07, 

where  K  is  the  block -diagonal  NlxNI  dimensional  matrix  defined  by 


K(O)  0  0 

0  /((  1 )  0 

0  0  0 


K{1-  1 


The  i-th  diagonal  block  K(i)  of  K  is  the  covariance  matrix  of  the  reflectance  process  at  the  i -th 
delay  bin.  The  imaging  problem  in  terms  of  these  expressions  is,  then,  to  estimate  the  matrix  K 
of  (27)  given  the  data  vector  r  of  (16).  The  scattering  function  then  can  be  determined  from  K 
by  using  (15). 


fi 


3.  Maximum- Likelihood  Imaging  for  the  Incomplete- Data  Model 

For  reasons  that  will  become  evident  in  the  next  section,  we  term  the  /V-dimensional  data 
vector  r  the  incomplete  data  for  the  radar- imaging  problem.  The  model  given  in  the  last  section 
for  this  incomplete  data  is  that  r  is  normally  distributed  with  zero  mean  and  covariance 
specified  in  (26).  Given  the  incomplete  data,  we  wish  to  estimate  the  covariance  K  of  the 
reflectance  process,  as  defined  in  (27).  To  do  this,  we  adopt  the  maximum- likelihood  method 
of  statistics,  which  selects  K  to  maximize  the  incomplete-data  toglikelihood 

LjK)~-Un{det[2ErS'  KS*N0l))-[-r'[2ErS‘  KS*  N0l)''  r. 

2  (28) 

where  the  maximization  is  subject  to  the  constraint  that  K  be  an  admissible  matrix,  where  by  an 
admissible  matrix  we  mean  a  matrix  having  the  block-diagonal  form  in  (27)  with  each  diagonal 
block  being  a  Hermitian-symmetric,  positive-semidefinite  Toeplitz- matrix. 

We  now  derive  a  necessary  condition,  termed  the  trace  condition ,  which  the  matrix 
maximizer  of  the  incomplete-data  toglikelihood  (28)  must  satisfy.  In  principle,  this 
equation  specifies  the  maximum-likelihood  estimate  of  K.  If  K  maximizes  L^(K)  of 
(28),  then  Lid(K+6K )  <r  Lid(K)  for  SK  small.  Equivalently,  the  first  derivative  is  zero, 

Lld[K*a6K)-Lld(K)  „  (29) 

,im-» - rsi - °' 

for  all  matrix  variations  SK  with  a  real  such  that  K  *  a SK  is  admissible.  As  shown  in  the 
Appendix,  this  implies  the  trace  condition, 

Tt[{2E  TS‘  KS  +  N  0l)~'{rr'  -  2E  TS'  KS  -  N  0l)(2ETS'  KS  +  N0l)~'  S'6Ks)-0,  (30) 

which  must  be  satisfied  by  the  maximum-likelihood  estimate  K.  Burg,  Luenberger,  and  Wenger 
[17)  have  studied  an  equivalent  problem  of  Toeplitz-constrained  covariance-estimation  and  have 
derived  the  trace  condition  using  a  different  approach. 

There  are  NI  unknowns  in  K.  Since  SK  must  be  a  block  diagonal  matrix  of 
Hermitian-symmetric  Toeplitz-matrices,  there  are  NI  parameters  in  SK  that  can  be  varied. 

These  variations  in  the  trace  condition  (30)  generate  NI  equations  in  the  unknown  elements  of 
K.  Thus,  in  principle,  the  trace  condition  produces  enough  equations  to  determine  the 


■v 

$ 
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maximum-likelihood  estimate  K.  However,  the  equations  are  complicated  due  to  the  inverse 
matrices  appearing  in  (30),  so  it  does  not  appear  to  be  feasible  to  determine  K  directly  from  the 
trace  condition,  which  motivates  our  development  of  the  iterative  approach  in  the  next  section. 
The  trace  condition  is  only  a  necessary  condition  which  the  estimate  K  must  satisfy.  For  it  to 
be  sufficient  as  well,  the  second  derivative  must  be  negative  along  all  admissible  variational 
directions  &K.  In  the  Appendix,  an  expression  for  the  second  derivative  in  the  direction  SK  is 
obtained. 

In  the  next  section,  we  will  develop  an  iterative  procedure  for  determining  a  sequence  of 
estimates  that  increase  the  likelihood  at  each  iteration  stage.  We  will  demonstrate  that  the  limit 
point  of  the  iterations  satisfies  the  trace  condition  (30). 


4.  Maximum- Likelihood  Imaging  for  the  Incomplete/Complete  Data  Model 

That  the  trace  condition  (30)  cannot  be  solved  directly  for  the  maximum- likelihood 
estimate  of  K  motivates  the  indirect  approach  we  now  take  of  embedding  the  imaging  problem 
in  a  larger,  seemingly  more  difficult  problem.  The  result  will  be  an  iterative  algorithm  which 
when  implemented  produces  a  sequence  of  admissible  matrices  K(°),  £(*),  ....  JCOO,  ...  having  the 
property  that  the  corresponding  sequence  of  incomplete-data  loglikelihoods  Z-idlATf0)],  L^A^1)], 

...,  £-id£AT(k)].  ...  is  nondecreasing  at  each  stage. 

Fuhrmann  and  Miller  [18]  have  recently  shown  that  maximum- likelihood  estimates  of 
Toeplitz-constrained  covariances  which  are  positive  definite  do  not  always  exist  when  given 
only  one  data  vector  r.  A  necessary  and  sufficient  condition  for  the  likelihood  function  to  be 
unbounded,  and  therefore  for  no  maximum-likelihood  estimate  to  exist,  is  that  there  be  a 
singular  Toeplitz  matrix  with  the  data  in  its  range  space.  For  our  imaging  problem,  this 
condition  is  that  there  exists  an  admissible  K  with 

2ETS'  KS  +  N0l 

singular  such  that 

r~[2ETS'  KS*N0l)a 

for  some  complex-valued  vector  a.  In  fact,  without  constraining  K  further  than  being  Toeplitz, 
a  sufficient  condition  that  a  singular  estimate  for  K  be  obtained  is  that  N0  -  0  and  there  exists  a 
singular  K  with  r  in  the  range  space  of  2E?S+KS.  The  argument  for  this  mirrors  that  of 
Fuhrmann  and  Miller  in  [18,  Theorem  1],  but  is  applied  to  the  complete  data  loglikelihood  (35). 
Furhmann  and  Miller  also  showed  that  even  if  the  true  covariance  had  eigenvalues  bounded 
from  above  and  below,  the  probability  that  there  exists  a  singular  Toeplitz  matrix  with  the  data 
in  its  range  can  be  very  close  to  one.  By  restricting  the  search  to  Toeplitz  matrices  with 
circulant  extensions,  they  were  able  to  show  that  the  probability  a  singular  circulant  Toeplitz 
matrix  has  the  data  in  its  range  space  is  zero.  Thus,  in  order  for  maximum-likelihood  estimates 
to  be  nonsingular  with  probability  one  for  all  nonnegative  values  of  No,  we  restrict  the  class  of 
admissible  Toeplitz  matrices  to  be  those  with  circulant  extensions  of  period  P,  where  P  is  equal 
to  or  greater  than  the  number  N  of  data  samples  available,  P  >  N.  What  we  envision  in 
adopting  this  constraint  is  that  for  each  delay  i,  the  N  sample  values  of  the  reflectance  b(n.i).  n 


*  0,  1 . N-I,  are  from  a  stationary  process  that  is  periodic  with  period  P ,  where  P  could  be 

some  large  but  finite  value.  These  N  sample  values  enter  the  incomplete  data  r  according  to 
(15)  and  (22).  By  using  the  expectation-maximization  algorithm  of  Dempster,  Laird,  and  Rubin 
[19],  we  shall  develop  a  sequence  of  admissible  matrices  that  have  the  maximum-likelihood 
estimate  of  K  subject  to  this  circulant  extension  as  a  limit  point.  The  approach  parallels  that  of 
Miller  and  Snyder  [15]  for  estimating  the  power  spectrum  of  a  time-series  from  a  single  set  of 
data.  An  important  benefit  of  introducing  the  periodic  extension  and  using  the  expecta¬ 
tion-maximization  algorithm  is  that  estimates  of  both  the  scattering  function  and  the  reflectance 
process  are  obtained  simultaneously  and  can  be  readily  viewed  as  target  images  in  range  and 
cross-range  coordinates;  thus,  the  procedure  proposed  may  be  considered  to  be  natural  for  the 
imaging  problem  because  both  types  of  images  considered  separately  in  the  past  are  obtained 
directly.  For  completeness,  we  also  include  in  the  Appendix  the  equations  obtained  using  the 
expectation-maximization  algorithm  for  estimating  general  Toeplitz  matrices  when  the 
assumption  of  a  circulant  extension  is  not  made. 

We  shall  introduce  a  modification  of  our  notation  to  indicate  that  the  N  samples  of  the 
reflectance  process  are  from  a  stationary  periodic-process  of  period  P.  Thus,  let  b^(i)  denote 
the  W-dimensional  vector  b(i)  of  (21).  We  now  think  of  *n(0  as  an  N-dimensional  subvector  of 
the  P-dimensional  vector  bp(i)  formed  from  samples  of  the  reflectance  process  over  a  full 
period, 

b  ( 0 ,  m  ♦  i ) 
b[  1 ,  m  *  t) 

bf  N  -  1  ,  m  ♦  i) 

6  f  P  -  1  ,  m  +  i ) 

If  /n  is  the  NxN  identity  matrix,  and  if  J  is  the  PxN  matrix  defined  by 


then  bn(i)  -  7+bp</).  Also,  let  b n  denote  the  A//-dimensional  vector  b  of  (20),  and  let  bp  be  the 
^/-dimensional  vector  with  i-th  block  element  bp(i).  Then,  bn  *  M+bp ,  where  M  is  the  PIxNI 
block  diagonal-matrix 


0  0 
J  0 


Vo  0  0 


(33) 


Furthermore,  let  Kn(i)  denote  the  NxN  Toeplitz  covariance- matrix  K(i)  of  by(i)  defined  in 
(25),  and  let  Kp(i)  denote  the  PxP  circulant  covariance- matrix  of  bp(i).  Then,  the  Toeplitz 
matrix  Kn(i)  is  the  upper  left  block  of  the  circulant  matrix  Kp(i),  as  given  by 

Kn  -  J'Kr[i)J. 


Lastly,  let  Kp  denote  the  PIxPl  block-diagonal  matrix  in  the  form  of  (27)  with  the  i-th  diagonal 
block  being  Kp(i).  Then,  if  Kn  denotes  the  NIxNI  matrix  K  of  (27),  there  holds 

Kn - 

For  use  with  the  expectation-maximization  algorithm,  we  identify  the  complete  data  as 
(6p,w),  where  w  is  the  W-dimensional  noise  vector  defined  in  (16).  We  note  from  (15),  (22), 
and  the  above  definitions  that  the  incomplete  data  r  can  be  obtained  from  the  complete  data 
according  to  the  mapping 


r  -  J2E  TS'  M'b,  +  ui.  (34) 

The  loglikelihood  function  Lc(j(/Cp)  of  the  complete  data  is  given  by 

LjKF)--^\n{de\.{Kr))-'-b',K-F'b,,  ('>5> 

where  all  terms  that  are  not  a  function  of  Kp  have  been  suppressed. 

Let  W  denote  the  PxP  discrete  Fourier-transform  matrix  scaled  so  that  the  columns  are 


orthonormal. 


(36) 
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where  wP  -  txp{-j2*/P).  Also,  let  lf>  be  the  PIxPl  block-diagonal  matrix 
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(37) 


Then,  frp  can  be  represented  in  rotated  coordinates  according  to 


-  1/ p b p  ' 


/  c  ( 0 ) 
'  c  ( 1 ) 


\c(/-iL 


(38) 


where  c(i)  »  lffep(/).  The  assumption  that  dp (/)  originates  from  a  periodic  process  implies  that 
the  /’/-dimensional  vector  crP  is  normally  distributed  with  zero  mean  and  diagonalized 
covariance 


E  p~l{cpc'p)-W  pKrW'p. 


(39) 


We  will  denote  the  (p+il)-//t  diagonal  element  of  £P  by  ov2(i)\  this  is  the  p-th  diagonal  element 
of  the  PxP  diagonal  matrix  E[c(j)c+(/)]. 


Substituting  the  expression 


Kp-WpZpWp 


(40) 


into  (35)  indicates  that  the  complete-data  loglikelihood  can  alternatively  be  expressed  as  a 
function  of  EP  according  to 


-  16 


(41) 


-~in(det(r,))-±c;r;'c, 

--Z  Z^lepWl-iZ  i>p(i')r*;2(*)- 

(-0  p-0  z  (- 0  p-0 


where  cp(i)  is  the  p-r/r  element  of  the  P-dimensional  vector  cp(»). 

The  expectation-maximization  algorithm  for  estimating  the  covariance  of  the  reflectance 
process  Kp  from  the  incomplete  data  r  is  an  alternating  maximization  procedure  in  which  a 
sequence  of  estimates  of  Ep  having  increasing  likelihood  is  obtained  first.  If  Ep(k)  denotes  the 
estimate  of  Ep  at  stage  k,  then  there  is  a  corresponding  element,  Kp(k)  -  fTp+EpOOlKp  in  a 
sequence  of  estimates  of  Kp  having  increasing  likelihood.  Likewise,  to  the  k-th  element  Ap(k)  of 
the  sequence  of  estimates  of  Kp,  there  is  a  corresponding  element,  AjfOO  m  Af+Ap(k)Af,  in  a 
sequence  of  estimates  of  having  increasing  likelihood. 

Each  iteration  stage  of  the  expectation-maximization  algorithm  has  an  expectation  (E)  step 
and  a  maximization  (M)  step  that  must  be  performed  to  get  to  the  next  step.  The  E-step 
requires  evaluation  of  the  conditional  expectation  of  the  complete-data  loglikelihood  (41)  given 
the  incomplete  data  r  and  assuming  that  the  covariance  defining  the  complete  data  is  EpOO, 

(42) 

From  (41),  we  have  that 

Z  in(ap(i))~Z 

p-0  ^  1-0 


r- 1 


Z  Etlcp(*)l 


\r,r[rk']o-p2{ 


p-0 


<-0 


The  M-step  yields  the  estimate  Ep(k+1)  at  stage  k+1  as  the  choice  of  Ep  that  maximizes  this 
conditional  expectation. 


ZTp*11  -  arg  max  [QCr^lr1,*1)], 

subject  to  the  constraint  that  the  maximizer  be  a  diagonal  covariance- matrix.  From  (43),  this 
maximization  yields  the  diagonal  matrix  Ep(k+1)  with  (p+iI)-/A  diagonal  element 


(a2(i))l**"-E[|cp(i)|2|r,r1;1]. 


(45) 


17  - 


Thus,  we  may  write  Ep(k+1)  as 


r'/.fc*,l-E[c,c;|r.r1/1], 


(46) 


where  the  "d"  over  the  equal  sign  means  that  the  diagonal  terms  in  the  matrix  on  the  left  side 
equal  the  diagonal  terms  in  the  matrix  on  the  right  side  and  that  all  the  off  diagonal  elements 
on  the  left  side  are  zero. 

The  above  expression  (46)  appears  to  be  complicated  because  of  the  several  matrices  we 
have  defined,  but  it  produces  a  sequence  of  covariance  estimates  having  a  straightforward 
interpretation.  If  we  form  the  matrix  Kp(k+t)  according  to 


then  we  find  that 


(47) 


K[fk'" 


0)  0  0 

0  1 )  0 

0  0  0 


0 

0 


(48) 


*!?•"(/- 1). 


where  Kp(k+1)(i)  is  a  PxP  circulant  matrix  interpreted  as  the  estimate  at  stage  k+1  of  the 
covariance  Kp(i)  of  the  /’-periodic  reflectance  process  at  delay  m+i.  Miller  and  Snyder  [15] 
show  that  the  (n,m)-element  of  this  circulant  matrix  is  given  by 


4  Z  E[t>(p,i)b*((p  +  n-m),.t)|r,/<',r*1], 

r  P-0 


(49) 


where  «z>p  ■  a  mod  P.  Equation  (49)  has  an  intuitively  appealing  form.  If  the  reflectivity 

process  b(n,i)  could  be  observed  for  all  instants  n  -  0,  1 . P-1  in  a  period  and  for  each  i 

independently,  then  the  maximum-likelihood  estimate  of  the  covariance  ATp(i)  would  be  the 
arithmetic  average  of  the  lagged  products 

(50) 


4  Z  t»(p.‘)t»*((p  +  n-m),,i). 
r  p.0 
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Equations  (48)  and  (49)  indicate  that  one  should  simply  substitute  the  conditional  mean  estimate 
of  an  unknown  lagged  product  into  this  expression  to  form  the  maximum- likelihood  estimate  of 
the  covariance  when  only  the  incomplete  data  are  known. 


* 


i 


■Sfl 

iyj 

t 


W 


estimating  Ep  and  Kp 

The  maximum-likelihood  estimate  of  Ep  is  a  limit  point  of  the  sequence  defined  in  (46). 
The  terms  on  the  right  side  of  this  equation  can  be  evaluated  as  follows.  Let  the 
conditional-mean  estimate  of  cp  in  terms  of  the  incomplete  data  r  be  defined  at  stage  k  by 

c,,*1-E[c/.|r,r‘.t|].  (5 

Then,  (46)  can  be  rewritten  in  the  form 

r1/* 1 1  -  E[(  c ,  -  c1;1 )  ( C.,  -  C )  ■ *  1  r ,  r  <;'  ] + e'/'e'/1* . 

(5: 

Now  estimating  cp  from  r,  where  from  (34)  and  (38) 

r-  J2E T S' M'W\c,  +  w,  (5; 

is  a  standard  problem  in  linear  estimation-theory.  From  Tretter  [20,  Ch.  14],  for  example,  we 
find  that 

elpl  -  r  */’  W  ,  MS  [  2  E  T  S *  M  *  w ;  r  l?] V  r  MS  -  N  0  I ] ' '  r . 

(5- 

Furthermore,  the  first  term  on  the  right  side  in  (32)  is  the  covariance  of  the  estimation  error 
when  cp  is  estimated  from  r.  Also  from  Tretter  [20,  Ch.  14],  we  have 

e[(c,-  cl;l)(c,-c,;,)*k.r1;1] 

-  r -  2 E  T  r W  , MS  [  2 E  T  S  *  M  *  W  ; r if1  w , MS  *  N  o /  ]  * 1  s  *  \t  *  w  ;r (,*' . 

In  summary,  the  following  steps  are  performed  to  produce  a  sequence  Ep(°),  EpU), 

...,  Ep(k),  ...  of  estimates  of  Ep  for  which  the  corresponding  sequence  of  likelihoods  is 
nondecreasing: 


r,  s-  ,  t  .  v  s  f 


s .  v  »  » •  >  a**  %  •  'i  9  v 


•  /  ;  ;  i1*,-  »'x  .*  «-<■'  v 


1 .  set  k  -  0,  select  a  starting  estimate  Ep(°); 

2.  calculate  the  estimate  of  cp  according  to  (54); 

3.  calculate  the  error  covariance  according  to  (55); 

4.  update  the  estimate  of  Ep  according  to  (52); 

5.  if  "last  iteration*  then  stop,  else  replace  k  by  k+1  and  go  to  2. 

The  starting  value  in  step  1  can  be  any  positive-definite,  diagonal  covariance- matrix  of 
dimension  PIxPI.  Clearly,  the  processing  indicated  in  (52)-(55)  is  a  nonlinear  function  of 
the  data. 


From  (40),  a  sequence  of  estimates  of  Kp  having  increasing  likelihood  is  obtained  from  the 
sequence  of  estimates  of  Ep  according  to  the  following  formula: 


K 


iti 

r 


(56) 


estimating  the  scattering  function 

Recall  that  /Cp  is  a  block-diagonal  matrix  with  the  i-th  diagonal  block  equal  to  the 
circulant  covariance  matrix  Kp(i)  of  the  P- periodic  reflectance  process  at  delay  m+i.  The  first 
column  of  this  matrix  is  given  by  Kp(i)e,  where  e  is  the  P-dimensional  unit  vector 


f'o\ 

0 

\0/ 


(57) 


Denote  the  scattering  function  of  the  P-periodic  reflectance  process  at  delay  m+i  by  S(i).  This 
is  a  P-dimensional  vector  with  p -th  element  given  from  (14)  by 


5(p-‘)  "  t)exp(-  ‘2nnp^ 

m  '[PWK r(i)e. 


j~ 


(58) 
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From  this  expression,  we  see  that  the  vector 


(59) 


is  a  vector  of  /  vectors  of  dimension  P,  the  i -th  of  which  is  the  scattering  function  at  delay  bin 
m+i.  Consequently  from  (56),  a  sequence  of  estimates  of  the  scattering  function  having 
increasing  likelihood  is  obtained  from  the  sequence  of  estimates  of  Dp  according  to  the  formula: 


fpr{!] t/J  • 

U 


(60) 


Now,  the  vector 


{FW 


e 

u 


appearing  in  (60)  is  an  /P-dimensional  vector  of  all  ones.  As  a  result,  the  quantity  in  (60)  is  a 
vector  whose  elements  are  the  diagonal  elements  of  Ep(k).  Thus,  at  iteration  stage  k,  the 
estimate  of  the  scattering  function  at  delay  m+i  is  given  by  the  P  diagonal  elements  of  the  i-th 
PxP  diagonal  block  of  Ep(k).  We  may,  therefore,  simply  regard  EpOO  as  the  stage  k  estimate  of 
the  scattering  function.  If  the  ( p+if)-th  diagonal  element  of  Ep(k)  is  the  placed  in  the  (p.i) 
element  of  a  Px/-dimensional  array,  as  p  varies  from  0  to  P-1  and  i  varies  from  0  to  1-1,  then 
the  result  may  be  displayed  as  the  target’s  scattering-function  image  at  stage  k  in  range  (/ 
coordinate)  and  cross  range  ( p  coordinate). 

i 
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estimating  the  reflectance  process 

It  is  interesting  to  note  that  the  k-th  stage  conditional-mean  estimate  of  cp,  given  the 
measurements  r  and  assuming  that  the  second-order  statistics  of  reflectance  are  given  by  the 
k-th  stage  estimate  of  the  scattering  function,  is  used  to  form  the  estimate  of  Spat  stage  k+1 
when  the  expectation-maximization  algorithm  is  used.  This  estimate  is  of  very  much  interest  in 
its  own  right  because,  from  (38),  cp (i)  is  the  Fourier  transform  of  the  reflectance  process  by(i). 
Thus,  if  the  (p+il)-th  element  of  this  estimate  is  placed  in  the  (p.i)  element  of  a  Pxl 
dimensional  array,  as  p  varies  from  0  to  P-1  and  »  varies  from  0  to  1-1 ,  then  the  result  may  be 
displayed  as  the  target’s  reflectance  image  at  stage  k  in  range  (»  coordinate)  and  cross  range  (p 
coordinate). 


convergence  issues 

There  are  some  important  properties  of  the  iteration  sequence  (46)  which  are  worth 
mentioning.  First,  each  step  is  in  an  improving  direction.  This  is  shown  by  writing  (32)  out  as 


ri*->  n 
A  r 


t-l*)  . 

**  r 


r'/'e 


,(*|  £-U)( 


where 


(61) 


6[k'  -2E  TW  ,MSK[?'y[rr'  -2E  TS'  fMS  -  N  0I)K{?'x  S'  M'W',., 


(62) 


and  where 

K[rki  -2F TS'  KTW,Zy?W fMS  ♦  N  0I 

is  the  k-/A  estimate  of  the  covariance  Kt  of  r.  Next,  the  trace  condition  (30)  which  the 
maximum- likelihood  estimate  must  satisfy  is  reexamined.  From  the  assumption  of  the 
P-periodicity  of  the  reflectance  process  and  the  matrix  definitions  given,  the  admissible 
variations  SK  must  be  of  the  form 


6K- M'W'r6£W  rM . 

Here,  $E  is  a  diagonal  matrix  of  the  same  dimensions  as  E.  The  trace  condition  (30)  then 
becomes 


Tt{k;'{2E  TS'  M'w;r  ,W  +  N0I -rr')K;'  S'  hrv;6EW 


(65) 


Using  the  fact  that  Tr(/4fl)  -  Tt{BA)  and  evaluating  this  trace  at  the  k-th  iterate,  we  see  that  the 
trace  on  the  left  side  of  (65)  is  equal  to 


(2£r)',Tr(e“l6r). 

According  to  (61),  Ep(k)  is  changed  at  each  stage  by  adding  the  diagonal  elements  of 


to  Ep(k).  Define 


r 


(*> 

r 


as  these  diagonal  elements.  Then,  evaluating  the  trace  at  this  variation  gives 


Tr^sr1*’)^. 

This  shows  that  the  variation 


(66) 


(67) 


(68) 


(69) 


sr141 

is  in  an  improving  direction.  Furthermore,  we  are  guaranteed  that  the  incomplete-data 
loglikelihood  is  nondecreasing  as  a  result  of  the  M-step  of  the  expectation-maximization 
algorithm.  At  this  step,  the  conditional  expectation  of  the  complete-data  loglikelihood  given  the 
incomplete  data  and  the  last  iterate  for  Ep  is  maximized  over  Ep.  As  shown  in  [15]  and  [19], 
this  implies  that  the  incomplete-data  loglikelihood  is  nondecreasing. 

In  the  Appendix,  we  show  that  if  N0  >  0  and  the  initial  guess  for  Ep  is  positive  definite, 
then  each  succeeding  guess  is  also  positive  definite.  This  gives  another  interpretation  of  (69) 
when  No  >  0.  Since  the  diagonal  elements  of  Ep(k)  are  positive,  (69)  holds  with  equality  if  and 
only  if  the  diagonal  elements  of  8(k)  are  zero.  It  follows  from  this  expression  (69)  that  a  second 
property  of  the  sequence  (46)  is  that  all  stable  points  of  the  sequence  satisfy  the  necessary  trace 
condition  (30).  This  is  easily  shown  by  noticing  that  if  Ep(k+i)  -  EpOO,  then  the  diagonal 
elements  (68)  are  zero.  But,  this  implies  that  the  diagonal  elements  of  6(k)  are  zero  and  hence 
that 


Tr(e|4|62:)-0 


(70) 
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for  all  diagonal  S£. 


computational  considerations 

The  computations  required  to  produce  radar  images  with  our  method  are  specified  by  (52), 
(54),  (55),  and  (60)  for  the  scattering-function  image  and  by  (52),  (54),  and  (55)  for  the 
reflectance  image.  The  number  of  iterations  of  these  equations  that  are  required  to  produce  an 
image  near  the  convergence  point  is  presently  unknown.  Our  experience  in  using  an  iterative 
algorithm  to  produce  maximum- likelihood  images  for  emission  tomography  suggests  that  50  to 
100  iterations  may  be  necessary,  but  this  is  only  a  guess  that  will  not  be  verified  until  some 
experiments  are  completed.  Some  form  of  specialized  processor  to  accomplish  each  iteration 
stage  efficiently  will  probably  be  needed  to  produce  images  in  practically  useful  times.  One 
possible  approach  is  the  following.  The  matrix  product 

r-j2ETW  fMS 

is  required  at  each  iteration  stage  and  does  not  change.  This  /YxJV-dimensional  matrix  can, 
therefore,  be  computed  once  off  line,  stored,  and  then  used  as  needed  during  on-line 
computations.  Then,  at  iteration  stage  k ,  the  following  on-line  computations  can  be  performed: 

1.  compute  the  NxN  matrix  A  defined  by  A  -  r+EP(k)r  +  Ntf; 

2.  compute  the  PIxN  matrix  B  defined  by  B  -  EpOOr; 

3.  compute  BA~lr  and  the  diagonal  elements  of  Ep(k)  -  BA-lB+. 

The  computations  in  3  can  be  accomplished  in  about  4N+PI-2  time  steps  using  the  systolic  array 
described  by  Comon  and  Robert  [21]  augmented,  as  they  suggest,  by  one  row  to  accomplish  the 
postmultiplication  of  ZM-i  by  r  and  by  IP  rows  to  accomplish  the  postmultiplication  by  B+.  The 
matrix  multiplications  in  1  and  2  for  determining  A  and  B  can  also  be  performed  rapidly  on  a 
systolic  array.  More  study  of  implementation  approaches  is  needed,  but  it  does  not  appear  that 
the  computational  complexity  of  our  new  imaging  algorithm  needs  to  be  a  limitation  to  its 
practical  use. 
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The  choice  of  N  and  7  is  important  for  the  computations.  These  parameters  are  selected  to 
achieve  a  desired  range  and  cross-range  resolution  and  are,  therefore,  problem  dependent,  but 
the  same  considerations  used  with  other  approaches  to  radar  imaging  can  be  used  in  selecting 
them.  On  the  other  hand,  the  choice  of  P  is  unique  to  our  approach.  As  stated,  we  need  P  > 
N,  but  no  upper  limit  is  given.  In  [18],  it  is  shown  that  as  P  increases  towards  infinity  so  does 
the  maximum  value  of  the  incomplete-data  loglikelihood  function,  with  probability  one.  Thus, 
P  cannot  be  made  arbitrarily  large  from  a  theoretical  standpoint.  Practically,  it  is  desired  to 
have  P  as  small  as  possible  to  decrease  the  memory  requirements  and  the  numerical  operations. 
The  natural  question,  then,  is  what  would  it  mean  to  have  P  equal  to  its  smallest  allowed  value 
For  this  choice,  IF  is  an  NxN  matrix  and  the  matrices  K( i )  are  Hermitian  symmetric, 
circulant  covariance  matrices.  Clearly,  for  a  process  which  is  not  circulant,  this  produces  an 
undesirable  estimate  for  K.  As  P  increases,  the  block  diagonal  matrices  in  KN  look  less  like 
circulant  matrices.  Hence,  there  is  a  tradeoff  involved  in  the  choice  of  P  between  the  practical 
constraints  of  storage  and  computations  and  desirable  estimates  for  K.  So  far  as  we  have  been 
able  to  determine,  there  is  no  theoretical  basis  for  selecting  P. 
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5.  Conclusions 


The  expressions  we  have  obtained  for  forming  images  of  diffuse,  fluctuating  radar  targets 
are  based  on  the  model  stated  in  Section  2.  The  target  reflectance  is  assumed  to  introduce 
wide-sense-stationary  uncorrelated-scattering  (WSSUS)  of  the  transmitted  signal  with  no  glint  or 
specular  components  being  present.  The  reflectance  process  is  assumed  to  be  a  WSSUS  Gaussian 
process  with  unknown  second-order  statistics  given  by  a  delay-dependent  covariance  or 
scattering  function.  Echos  of  the  transmitted  signal  are  received  from  all  the  reflecting  patches 
that  make  up  the  target.  Each  patch  introduces  some  propagation  delay,  doppler  shift,  and 
random  amplitude-scaling  into  the  signal  it  reflects.  The  superposition  of  the  echos  from  all  the 
patches  is  received  in  additive  noise.  Thus,  the  reflectance  process  is  only  observed  indirectly 
following  a  linear  superposition  and  in  additive  noise,  so  neither  the  reflectance  process  nor  its 
second-order  statistics  are  known.  Target  images  are  made  by  displaying  estimates  of  either  the 
reflectance  process  or  its  second-order  st?fistics  (scattering  function)  based  on  processing  the 
received  signal.  In  Section  3,  we  derived  the  trace  condition  which  the  maximum-likelihood 
estimate  of  the  covariance  of  the  reflectance  must  satisfy,  and  we  concluded  that  this  condition 
is  too  complicated  to  solve  explicitly  for  the  estimate.  This  motivated  the  introduction  in 
Section  4  of  the  incomplete-complete  data  model  and  the  use  of  the  expectation-maximization 
algorithm,  which  results  in  a  sequence  of  estimates  of  the  scattering  function  having  increasing 
likelihood.  A  corresponding  sequence  of  estimates  of  the  reflectance  process  is  also  obtained. 

There  are  a  number  of  issues  yet  to  be  resolved  for  the  approach  we  have  presented,  and 
we  are  addressing  these.  Glint  and  specular  components  in  the  return  echos  need  to  be 
accommodated.  The  selection  of  transmitted  signals  to  produce  good  images  is  an  important 
subject  about  which  little  study  has  been  made.  The  quality  of  target  images  obtained  with  our 
new  approach  is  not  known  at  present;  to  study  this  issue,  we  are  presently  implementing  a 
computer  simulation  so  that  comparisons  to  alternative  processing  strategies  can  be  made.  The 
equations  we  have  developed  are  computationally  demanding,  so  special  processing  architectures 
will  be  important  to  make  their  use  practical. 


6.  References 


1.  M.  J.  Prickett  and  C.  C.  Chen,  "Principles  of  Inverse  Synthetic  Aperture  Radar  (ISAR) 
Imaging,”  IEEE  EASCON  Record,  pp.  340-345,  September  1980. 

2.  K..  Ruttenberg  and  L.  Chanzit,  "High  Range  Resolution  by  Means  of  Pulse  to  Pulse 
Frequency  Shifting,"  IEEE  EASCON  Record,  1968. 

3.  C.-C.  Chen  and  H.  C.  Andrews,  "Multifrequency  Imaging  of  Radar  Turntable  Data,"  IEEE 
Trans,  on  Aerospace  and  Electronic  Systems,  Vol.  AES-16,  pp.  15-22,  January  1980. 

4.  D.  L.  Mensa,  High-Resolution  Radar  Imaging,  Artech  House,  Dedham,  Mass.,  1984. 

5.  J.  L.  Walker,  "Range-Doppler  Imaging  of  Rotating  Objects,"  IEEE  Transactions  on 
Aerospace  and  Electronic  Systems,  Vol.  AES- 16,  No.  1,  pp.  23-52,  January  1980. 

6.  D.  C.  Munson,  Jr.,  J.  D.  O'Brian,  and  W.  K.  Jenkins,  "A  Tomographic  Formulation  of 
Spotlight- Mode  Synthetic  Aperture  Radar,"  Proceedings  of  the  IEEE,  Vol.  71,  pp.  917-925, 
August  1983. 

7.  R.  E.  Blahut,  "Segmented-Chirp-Waveform  Implemented  Radar-System,"  U.  S.  Patent  No. 
4309703,  January  1982. 

8.  M.  Bernfeld,  "CHIRP  Doppler  Radar,"  Proceedings  of  the  IEEE,  Vol.  72,  pp.  540-541, 
April  1984. 

9.  D.  L.  Snyder,  H.  J.  Whitehouse,  J.  T.  Wohlschlaeger,  and  R.  Lewis,  "A  New  Approach  to 
Radar/Sonar  Imaging,"  Proc.  SPIE  Conference  on  Advanced  Algorithms  and  Architectures, 
Vol.  696,  San  Diego,  CA,  August  1986. 

10.  D.  L.  Snyder,  L.  J.  Thomas,  Jr.,  and  M.  M.  Ter-Pogossian,  "A  Mathematical  Model  for 
Positron  Emission  Tomography  Systems  Having  Time-of-Flight  Measurements,"  IEEE 
Transactions  on  Nuclear  Science,  Vol.  NS-28,  pp.  3575-3583,  June  1981. 

11.  D.  L.  Snyder  and  D.  G.  Politte,  "Image  Reconstruction  from  List-Mode  Data  in  an 
Emission  Tomography  System  Having  Time-of-Flight  Measurements,"  IEEE  Transactions 
on  Nuclear  Science,  Vol.  NS-30,  pp.  1843-1849,  1983. 

12.  H.  L.  Van  Trees,  Estimation,  Detection,  and  Modulation  Theory,  Vol.  3,  John  Wiley  and 
Sons. 

13.  J.  Shapiro,  B.  A.  Capron,  and  R.  C.  Harney,  "Imaging  and  Target  Detection  with  a 
Hetrodyne-Reception  Ootical  Radar,  Applied  Optics,  Vol.  20,  No.  19,  pp.  3292-3313, 
October  1981. 

14.  V.  S.  Frost,  J.  A.  Stiles,  K.  S.  Shanmugan,  and  J.  C.  Holtzman,  "A  Model  for  Radar  Images 
and  Its  Application  to  Adaptive  Digital  Filtering  of  Multiplicative  Noise,"  IEEE 
Transactions  on  Pattern  Analysis  and  Machine  Intelligence,  Vol.  PAMI-4,  No.  2,  pp. 
643-652,  March  1982. 

15.  M.  I.  Miller  and  D.  L.  Snyder,  "The  Role  of  Likelihood  and  Entropy  in  Incomplete-Data 
Problems:  Applications  to  Estimating  Point-Process  Intensities  and  Toeplitz-Constrained 
Covariances,"  Proceedings  of  the  IEEE,  Vol.  75,  No.  7,  July  1987. 

16.  M.  Turmon  and  M.  I.  Miller,  "Performance  Evaluation  of  Maximum-Likelihood  Estimates 
of  Toeplitz  Convariances  Generated  with  the  Expectation-Maximization  Algorithm,"  1987 
Conference  on  Information  Sciences  and  Systems,  Johns  Hopkins  University,  Baltimore, 
MD,  pp.  25-27,  1987. 


17.  J.  P.  Burg,  D.  G.  Luenberger,  and  D.  L.  Wenger,  "Estimation  of  Structured  Covariance 
Matrices,"  Proc.  IEEE,  Vol.  70,  No.  9,  pp.  963-974,  September  1982. 

18.  D.  R.  Fuhrmann  and  M.  I.  Miller,  "On  the  Existence  of  Positive  Definite  Maximum-Likeli¬ 
hood  Estimates  of  Structured  Covariance  Matrices,"  IEEE  Transactions  on  Information 
Theory,  in  review. 

19.  A.  D.  Dempster,  N.  M.  Laird,  and  D.  B.  Rubin,  "Maximum  Likelihood  from  Incomplete 
Data  Via  the  EM  Algorithm,"  J.  of  the  Royal  Statistical  Society,  Vol.  B,  39,  pp.  1-37,  1977. 

20.  S.  A.  Tretter,  Introduction  to  Discrete-Time  Signal  Processing ,  John  Wiley  and  Sons,  New 
York,  1976. 

21  P.  Comon  and  Y.  Robert,  "A  Systolic  Array  for  Computing  BA-1,"  IEEE  Trans,  on 

Acoustics,  Speech  and  Signal  Processing,  Vol.  ASSP-35,  No.  6,  pp.  717-723,  June  1987. 


-  28  - 


7.  Appendix 

derivation  of  the  trace  condition  ( 30) 

From  the  definition  of  the  loglikelihood  function  in  (28),  we  have 

“  (A 

"~^ar '{{Kr  +  a2ETS'6KS) ''-K'riy-~ln(det(Kr  +  a2ETS*6KS)dat(K;')), 

where  Kr  is  the  covariance  of  the  incomplete  data  r  as  given  in  (26).  Examining  the  first  term 
on  the  right,  we  have  that  it  equals 


I  +  a2ETS'6KSK'r'Yl- l)r 


K;'(a2E  TS'  6KSK'rl)(l  +  a2E  TS*  6KSK~rl)  ' r 

2a 


^r'K;'2ETS^6KSK;'r  +  0{a). 


Examining  the  second  term  on  the  right  in  (Al),  we  have 

-^-ln{det(l  +  a2FTS'6KS/Cr'}}  - -~ln{det{/ +aB)) 

2  a  2a 

-  -  “Inf  1  +  aTr(B)  +  ...  +  a*det(fi)) 

2a 

i 

--~Tr(B)  +  0(a), 

where  B  is  defined  in  the  first  equality.  Equations  (A2)  and  (A3)  imply  that 
r'  K;'2ETS'6KSK~r'r-Tr{2ETS'6KSK;')-0, 

which  is  the  trace  condition  (30). 

derivation  of  the  sign-definiteness  condition 

To  check  the  sign-definiteness  of  the  second  derivative  of  the  loglikelihood  L-^(K),  we 
form  the  limit 


Tr({2ETS'{K  +  a6K)S  +  N0l)'  [rr' -  N 0I -2E rS'{K  + a6K)s) 

x(2£  TS'KS  +  a2ETS'6KS  +  N0iylS*6KS 
-K'rl(rr'  -  N  0I -2E  TS'  KS)K~r'  S'  6KS) 
mTr{K'rl2E  TS'  6KSK~rl[2E  TS'  KS  +  N 0I  -2rr')K~l2E TS' 6 KS). 

(A5) 

A  necessary  condition  for  K  to  be  a  relative  maximum  is  that  this  last  expression  evaluated  at  K 
be  equal  to  or  less  than  zero  for  all  admissible  variations  SK.  Under  the  assumptions  in  Section 
4,  admissible  variations  are  given  by  (64).  Substituting  (64)  into  (A5)  and  evaluating  for  all' 
diagonal  matrices  5E  gives  the  second-order  necessary  condition.  A  sufficient  condition  for  K 
to  be  a  strict  local  maximum  is  that  the  trace  condition  (30)  is  satisfied  and  that  (A5)  is  strictly 
negative  for  all  admissible  variations. 

estimating  a  general  Toeplitz  matrix 

In  Section  4,  we  derived  a  sequence  of  estimates  for  a  covariance  matrix  subject  to  the 
constraint  that  the  estimates  must  be  circulant  Toeplitz  matrices.  For  completeness,  we  develop 
and  discuss  in  this  section  the  equations  for  estimating  a  covariance  matrix  subject  to  the 
weaker  constraint  that  the  estimates  be  general  Toeplitz  matrices.  Similar  equations  for  other 
constraints  on  the  Toeplitz  matrices  are  easily  obtained  by  mimicking  the  steps  in  the  main  body 
of  this  paper. 

Let  the  complete  data  be  {b,w),  and  let  b  be  normally  distributed  with  zero  mean  and 
covariance  K,  as  given  in  (27).  The  complete-data  loglikelihood  is  given  in  (35).  Maximizing 
this  function  gives  the  trace  condition 

Tr(K',(bb*-K)K'16A')-0,  (A6) 

which  the  maximum-likelihood  estimate  K  must  satisfy.  Performing  the  E  and  M  steps  of  the 
expectation-maximization  algorithm  yields  the  following  iteration  sequence  for  the  elements 
K(n,i).  n  •  0,  l . N-l,  of  the  covariance  matrix  K(i)  defined  in  (25): 

i  w-«-i  (A7) 

AT 1  ’( n , i)  —  — - E[  £  b{j,m  +  i)b*[j  +  n,m  +  i)\r,Kik)]. 


In  matrix  form, 
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l**i) . 
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(*i . 


2  EtK 


lk]SKik]'l[rr'  -2ETS'  K[k)S-  N0l)Klk]~'  S'  K 


iti 


where 


(A8) 


K[k]  ~2ETS'  K[k]S  +  N0I .  (A9> 

If  this  iteration  converges  to  a  stable  point,  then  the  trace  condition  is  satisfied  at  that  point,  as 
may  be  shown  by  using  the  same  arguments  as  in  Section  4.  It  is  worth  restating  that  the  reason 
this  iteration  is  not  recommended  here  is  that  the  probability  that  the  iteration  sequence 
generates  a  singular  estimate  for  K  approaches  one  as  N  gets  large.  By  restricting  consideration 
to  Toeplitz  matrices  with  circulant  extensions,  the  loglikelihood  function  is  bounded  with 
probability  one  for  finite  extensions  and  a  positive  definite  K  is  generated  with  probability  one, 
as  proven  by  Fuhrmann  and  Miller  [18]. 


proof  that  Ep(k)  is  positive  definite  for  every  k 

Assume  that  the  initial  guess  Ep(<0  for  £p  is  positive  definite  and  that  N0  >  0.  We  will  now 
show  that  if  EpOO  is  positive  definite,  then  so  is  Ep(k+1),  and  thus,  by  induction,  £p(k)  is 
positive  definite  for  all  k.  One  key  to  following  this  derivation  is  the  matrix  identity 

B[I  +  AB)~' -(I  +  BA)~'B.  <A1°) 

This  identity  is  used  to  rewrite  (62)  as 

E[P"]  -  H(N02EtL{;]W pMSS'  M'W'pE[p]*  NlL{;] 

+  2ETL[!]W  ,MSrr'  S' 

-  N 0 2  E T(  H  Elk 1 W F MS )  ( S '  M '  W \ L^k]  H ' ) 

+  2Et{hE[?W  FMSr){r'S'  M'W'PL[?H') 

+  NlHE[?H', 
where  we  have  defined  H  according  to 

H  -  ( 2  E  T  E  {Fl 1 W  F  U  S  S '  M '  W  F  +  N  0  / )  ‘ 1 . 


(All) 
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Clearly,  all  the  diagonal  elements  of  (All)  are  greater  than  or  equal  to  zero.  To  show  that  they 
are  strictly  positive,  we  look  at  the  last  term  and  get  that  the  i-th  diagonal  element  is 


IF-  1 


[nIhl"h')u  -n*  I  (H)M(rJf')  (//*) 

y-o 


y-o 


(A13) 


which  is  clearly  positive  when  So  >  0  since  H  is  invertible  and  all  diagonal  elements  of  Ep(k) 
are  positive. 


